From 39d62f5e99306fd709f15429df16309c1799d3d5 Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sat, 11 Jul 2026 19:19:46 +0200 Subject: [PATCH 01/17] Consolidate git submodules under Submodules/; add nanoflann, picoflann, fcpw Toward a benchmark suite (issue #109), gather all submodules under a single Submodules/ directory: - Move common-3d-test-models -> Submodules/common-3d-test-models, updating every relative mesh path in the Examples (MeshSDF, CSGUnion) and Integrations (AMReX MeshSDF/PaintEB) main.cpp defaults, Examples/CMakeLists.txt, and the docs (ObtainingEBGeometry/ExampleMeshSDF/ExampleCSGUnion.rst, READMEs) + the REUSE.toml note. URLs and the submodule's own name are left untouched. Verified MeshSDF and CSGUnion still load the mesh at the new path. - Add the benchmark comparison libraries as submodules: Submodules/nanoflann, Submodules/picoflann, Submodules/fcpw (all shallow). fcpw's Eigen dependency is its own nested submodule (deps/eigen). reuse lint compliant (submodule contents carry their own licensing). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01HS5LFQihKoMrE6nbPaQW9R --- .gitmodules | 14 +++++++++++++- Docs/Sphinx/source/ExampleCSGUnion.rst | 2 +- Docs/Sphinx/source/ExampleMeshSDF.rst | 2 +- Docs/Sphinx/source/ObtainingEBGeometry.rst | 6 +++--- Examples/CMakeLists.txt | 2 +- Examples/CSGUnion/README.md | 2 +- Examples/CSGUnion/main.cpp | 2 +- Examples/MeshSDF/README.md | 2 +- Examples/MeshSDF/main.cpp | 4 ++-- Examples/README.md | 2 +- Integrations/AMReX/MeshSDF/main.cpp | 2 +- Integrations/AMReX/PaintEB/main.cpp | 2 +- REUSE.toml | 2 +- .../common-3d-test-models | 0 Submodules/fcpw | 1 + Submodules/nanoflann | 1 + Submodules/picoflann | 1 + 17 files changed, 31 insertions(+), 16 deletions(-) rename common-3d-test-models => Submodules/common-3d-test-models (100%) create mode 160000 Submodules/fcpw create mode 160000 Submodules/nanoflann create mode 160000 Submodules/picoflann diff --git a/.gitmodules b/.gitmodules index 1915074c..ed5eb559 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,4 +1,16 @@ [submodule "common-3d-test-models"] - path = common-3d-test-models + path = Submodules/common-3d-test-models url = https://github.com/alecjacobson/common-3d-test-models shallow = true +[submodule "Submodules/nanoflann"] + path = Submodules/nanoflann + url = https://github.com/jlblancoc/nanoflann.git + shallow = true +[submodule "Submodules/picoflann"] + path = Submodules/picoflann + url = https://github.com/rmsalinas/picoflann.git + shallow = true +[submodule "Submodules/fcpw"] + path = Submodules/fcpw + url = https://github.com/rohan-sawhney/fcpw.git + shallow = true diff --git a/Docs/Sphinx/source/ExampleCSGUnion.rst b/Docs/Sphinx/source/ExampleCSGUnion.rst index 483e589b..12ab8cc4 100644 --- a/Docs/Sphinx/source/ExampleCSGUnion.rst +++ b/Docs/Sphinx/source/ExampleCSGUnion.rst @@ -16,4 +16,4 @@ for how to compile it with CMake, GNU Make, or a direct compiler invocation. cd Examples/CSGUnion ./CSGUnion.ex # defaults to cow.obj - ./CSGUnion.ex ../../common-3d-test-models/data/cow.obj + ./CSGUnion.ex ../../Submodules/common-3d-test-models/data/cow.obj diff --git a/Docs/Sphinx/source/ExampleMeshSDF.rst b/Docs/Sphinx/source/ExampleMeshSDF.rst index 59fb5af4..1f0863df 100644 --- a/Docs/Sphinx/source/ExampleMeshSDF.rst +++ b/Docs/Sphinx/source/ExampleMeshSDF.rst @@ -23,7 +23,7 @@ for how to compile it with CMake, GNU Make, or a direct compiler invocation. cd Examples/MeshSDF ./MeshSDF.ex # defaults to armadillo.obj - ./MeshSDF.ex ../../common-3d-test-models/data/cow.obj # or pick another mesh + ./MeshSDF.ex ../../Submodules/common-3d-test-models/data/cow.obj # or pick another mesh With no argument the example loads ``armadillo.obj`` from the ``common-3d-test-models`` submodule, so make sure it is checked out first (see :ref:`Sec:Cloning`). diff --git a/Docs/Sphinx/source/ObtainingEBGeometry.rst b/Docs/Sphinx/source/ObtainingEBGeometry.rst index 10217425..342960b5 100644 --- a/Docs/Sphinx/source/ObtainingEBGeometry.rst +++ b/Docs/Sphinx/source/ObtainingEBGeometry.rst @@ -13,7 +13,7 @@ Clone the repository from `GitHub `_: The core library is header-only and completely self-contained once cloned. However, the ready-to-run examples in :file:`Examples/` read surface meshes from the `common-3d-test-models `_ collection, -which is bundled as a git submodule (``common-3d-test-models/``) at the repository root. +which is bundled as a git submodule (``Submodules/common-3d-test-models/``) at the repository root. If you intend to run the bundled examples, clone with the submodule in one step instead: .. code-block:: bash @@ -26,13 +26,13 @@ If you already cloned without ``--recurse-submodules``, fetch the submodule afte git submodule update --init --recursive -The meshes are then available as ``.obj`` files under ``common-3d-test-models/data/``. +The meshes are then available as ``.obj`` files under ``Submodules/common-3d-test-models/data/``. Some mesh-based examples take a mesh path on the command line, resolved relative to the run directory (each example is run from its own source folder), for example: .. code-block:: bash - ./a.out ../../common-3d-test-models/data/armadillo.obj + ./a.out ../../Submodules/common-3d-test-models/data/armadillo.obj Running an example with no argument falls back to a default mesh from the submodule, so the submodule must be checked out for the examples to run. diff --git a/Examples/CMakeLists.txt b/Examples/CMakeLists.txt index 6defb562..86549f74 100644 --- a/Examples/CMakeLists.txt +++ b/Examples/CMakeLists.txt @@ -5,7 +5,7 @@ # and layer on the repo-side concerns: optional sanitizers and test registration. # # The binary is run with its source directory as the working directory so that -# relative mesh paths (e.g. "../../common-3d-test-models/data/*.obj") +# relative mesh paths (e.g. "../../Submodules/common-3d-test-models/data/*.obj") # in each example resolve correctly. function(ebgeometry_add_example name) add_subdirectory(${name}) diff --git a/Examples/CSGUnion/README.md b/Examples/CSGUnion/README.md index c5bcc62f..05ea1262 100644 --- a/Examples/CSGUnion/README.md +++ b/Examples/CSGUnion/README.md @@ -63,7 +63,7 @@ Running Run from this directory so the default mesh path resolves: ./CSGUnion.ex - ./CSGUnion.ex ../../common-3d-test-models/data/cow.obj + ./CSGUnion.ex ../../Submodules/common-3d-test-models/data/cow.obj With no argument the example loads `cow.obj` from the `common-3d-test-models` submodule, so make sure it is checked out diff --git a/Examples/CSGUnion/main.cpp b/Examples/CSGUnion/main.cpp index e7590371..15c299db 100644 --- a/Examples/CSGUnion/main.cpp +++ b/Examples/CSGUnion/main.cpp @@ -40,7 +40,7 @@ main(int argc, char* argv[]) // Mesh to merge with the sphere. Pass a path on the command line, or fall back to an OBJ file from // the common-3d-test-models submodule (path is relative to this example's source folder, where the // executable is run). See the "Building and using" docs for how to fetch the submodule. - std::string file = "../../common-3d-test-models/data/cow.obj"; + std::string file = "../../Submodules/common-3d-test-models/data/cow.obj"; if (argc >= 2) { file = std::string(argv[1]); } diff --git a/Examples/MeshSDF/README.md b/Examples/MeshSDF/README.md index c81d9b40..c60123ab 100644 --- a/Examples/MeshSDF/README.md +++ b/Examples/MeshSDF/README.md @@ -68,7 +68,7 @@ Running Run from this directory so the default mesh path resolves: ./MeshSDF.ex - ./MeshSDF.ex ../../common-3d-test-models/data/cow.obj + ./MeshSDF.ex ../../Submodules/common-3d-test-models/data/cow.obj With no argument the example loads `armadillo.obj` from the `common-3d-test-models` submodule, so make sure it is checked out diff --git a/Examples/MeshSDF/main.cpp b/Examples/MeshSDF/main.cpp index 874497b8..1380c2c4 100644 --- a/Examples/MeshSDF/main.cpp +++ b/Examples/MeshSDF/main.cpp @@ -31,9 +31,9 @@ int main(int argc, char* argv[]) { // Path to a surface mesh (STL/PLY/VTK/OBJ). Pass one on the command line, e.g. - // ./a.out ../../common-3d-test-models/data/cow.obj + // ./a.out ../../Submodules/common-3d-test-models/data/cow.obj // Paths are resolved relative to the run directory (this example's source folder when run via ctest). - std::string file = "../../common-3d-test-models/data/armadillo.obj"; + std::string file = "../../Submodules/common-3d-test-models/data/armadillo.obj"; if (argc == 2) { file = std::string(argv[1]); diff --git a/Examples/README.md b/Examples/README.md index bc454aa6..e7903cfd 100644 --- a/Examples/README.md +++ b/Examples/README.md @@ -17,5 +17,5 @@ Examples that couple EBGeometry to a third-party application code (AMReX, Chombo [`Integrations/`](../Integrations/README.md) at the repository root instead. The mesh files used by the examples come from the [common-3d-test-models](https://github.com/alecjacobson/common-3d-test-models) -git submodule (`common-3d-test-models/`) at the repository root. Clone the repository with `--recurse-submodules` (or run +git submodule (`Submodules/common-3d-test-models/`) at the repository root. Clone the repository with `--recurse-submodules` (or run `git submodule update --init --recursive`); see the "Building and using" documentation for details. diff --git a/Integrations/AMReX/MeshSDF/main.cpp b/Integrations/AMReX/MeshSDF/main.cpp index a46146d5..46a89b50 100644 --- a/Integrations/AMReX/MeshSDF/main.cpp +++ b/Integrations/AMReX/MeshSDF/main.cpp @@ -96,7 +96,7 @@ main(int argc, char* argv[]) // Mesh file (STL/PLY/VTK/OBJ). Override with 'filename=' in the inputs file. The default is an OBJ // from the common-3d-test-models submodule (see the "Building and using" docs for how to fetch it); the // path is relative to this example's source folder, where the executable is run. - filename = "../../common-3d-test-models/data/armadillo.obj"; + filename = "../../Submodules/common-3d-test-models/data/armadillo.obj"; pp.query("filename", filename); RealBox rb({-100, -75, -100}, {100, 125, 100}); // bounds the default armadillo mesh diff --git a/Integrations/AMReX/PaintEB/main.cpp b/Integrations/AMReX/PaintEB/main.cpp index bf7cee59..7cf1cf86 100644 --- a/Integrations/AMReX/PaintEB/main.cpp +++ b/Integrations/AMReX/PaintEB/main.cpp @@ -112,7 +112,7 @@ main(int argc, char* argv[]) // Mesh file (STL/PLY/VTK/OBJ). Override with 'filename=' in the inputs file. The default is a // small OBJ from the common-3d-test-models submodule (see the "Building and using" docs for how to // fetch it); the path is relative to this example's source folder, where the executable is run. - filename = "../../common-3d-test-models/data/suzanne.obj"; + filename = "../../Submodules/common-3d-test-models/data/suzanne.obj"; pp.query("filename", filename); RealBox rb({-5, -1, 2}, {0, 4, 6}); // bounds the default suzanne mesh diff --git a/REUSE.toml b/REUSE.toml index 99f62e58..9f53e573 100644 --- a/REUSE.toml +++ b/REUSE.toml @@ -6,7 +6,7 @@ # Source code carries inline SPDX headers instead; see Source/ and Examples/*/main.cpp. # # NOTE: third-party example meshes are no longer redistributed here; they are pulled in -# via the common-3d-test-models git submodule at the repository root (a separate repository +# via the common-3d-test-models git submodule under Submodules/ (a separate repository # with its own licensing, which REUSE does not descend into). version = 1 diff --git a/common-3d-test-models b/Submodules/common-3d-test-models similarity index 100% rename from common-3d-test-models rename to Submodules/common-3d-test-models diff --git a/Submodules/fcpw b/Submodules/fcpw new file mode 160000 index 00000000..61814ff0 --- /dev/null +++ b/Submodules/fcpw @@ -0,0 +1 @@ +Subproject commit 61814ff0c7b69d61dac3b9725cda1541b1b3ec4f diff --git a/Submodules/nanoflann b/Submodules/nanoflann new file mode 160000 index 00000000..ff0eb50a --- /dev/null +++ b/Submodules/nanoflann @@ -0,0 +1 @@ +Subproject commit ff0eb50a4b972c642c54fe453095b724d9ae7c60 diff --git a/Submodules/picoflann b/Submodules/picoflann new file mode 160000 index 00000000..d5fd165a --- /dev/null +++ b/Submodules/picoflann @@ -0,0 +1 @@ +Subproject commit d5fd165aa8f0f91ceb6277252db3e8f4455cb470 From d1124ff81a1881f740dd00728dd194c85448f84f Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sat, 11 Jul 2026 19:29:14 +0200 Subject: [PATCH 02/17] Add Benchmark/ suite: PointCloudBVH vs picoflann/nanoflann, TriMeshSDF vs fcpw Per issue #109, add a benchmark suite comparing EBGeometry against other geometry-query libraries (pinned as submodules under Submodules/). Illustrative and not CI-tested, like Integrations/. - Benchmark/NearestNeighbor: all-nearest-neighbor over a 500k point cloud -- PointCloudBVH vs picoflann vs nanoflann. Flann queries iterated in Hilbert order for a fair (warm-cache) comparison; the one-time sort they need is reported separately (EBGeometry reuses its build order for free). PointCloudBVH ~0.29 us/pt, nanoflann ~0.25, picoflann ~0.47; folding the sort in, PointCloudBVH wins the end-to-end job. - Benchmark/MeshSDF: closest-point on a triangle mesh (armadillo, ~100k tris, float) -- TriMeshSDF vs fcpw, the task both are built for. Mesh parsed once; each builds its own BVH over the same triangles and answers the same queries; distances cross-checked (0 mismatches). TriMeshSDF ~2.6 us/query vs fcpw ~7.5 -- with the caveat, documented, that fcpw runs its scalar Eigen fallback here (Enoki off) while TriMeshSDF is SIMD, so fcpw's vectorized path would narrow the query gap. Adds a Benchmark.rst docs page (wired into the toctree) and a REUSE.toml block for the benchmark build files/READMEs (main.cpp carry inline SPDX). Both benchmarks build via their GNUmakefiles against the Submodules/ deps. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01HS5LFQihKoMrE6nbPaQW9R --- Benchmark/MeshSDF/GNUmakefile | 25 +++ Benchmark/MeshSDF/main.cpp | 131 +++++++++++++++ Benchmark/NearestNeighbor/GNUmakefile | 24 +++ Benchmark/NearestNeighbor/main.cpp | 222 ++++++++++++++++++++++++++ Benchmark/README.md | 67 ++++++++ Docs/Sphinx/source/Benchmark.rst | 29 ++++ Docs/Sphinx/source/index.rst | 1 + REUSE.toml | 10 ++ 8 files changed, 509 insertions(+) create mode 100644 Benchmark/MeshSDF/GNUmakefile create mode 100644 Benchmark/MeshSDF/main.cpp create mode 100644 Benchmark/NearestNeighbor/GNUmakefile create mode 100644 Benchmark/NearestNeighbor/main.cpp create mode 100644 Benchmark/README.md create mode 100644 Docs/Sphinx/source/Benchmark.rst diff --git a/Benchmark/MeshSDF/GNUmakefile b/Benchmark/MeshSDF/GNUmakefile new file mode 100644 index 00000000..8ec70447 --- /dev/null +++ b/Benchmark/MeshSDF/GNUmakefile @@ -0,0 +1,25 @@ +# Benchmark: EBGeometry TriMeshSDF vs fcpw (closest-point on a triangle mesh). +# Requires the Submodules/fcpw submodule plus its Eigen dependency, and the mesh submodule: +# git submodule update --init Submodules/fcpw Submodules/common-3d-test-models +# git -C Submodules/fcpw submodule update --init deps/eigen # Eigen only (skip GPU slang-rhi) +# EBGEOMETRY_HOME defaults to the repo root (two levels up). + +EBGEOMETRY_HOME ?= ../.. +CXX ?= g++ +CXXFLAGS ?= -std=c++17 -O3 -march=native + +INCLUDES := -I$(EBGEOMETRY_HOME) \ + -I$(EBGEOMETRY_HOME)/Submodules/fcpw/include \ + -I$(EBGEOMETRY_HOME)/Submodules/fcpw/deps/eigen + +TARGET := MeshSDF.ex + +$(TARGET): main.cpp + $(CXX) $(CXXFLAGS) $(INCLUDES) $< -o $@ + +.PHONY: run clean +run: $(TARGET) + ./$(TARGET) + +clean: + $(RM) $(TARGET) diff --git a/Benchmark/MeshSDF/main.cpp b/Benchmark/MeshSDF/main.cpp new file mode 100644 index 00000000..02713739 --- /dev/null +++ b/Benchmark/MeshSDF/main.cpp @@ -0,0 +1,131 @@ +// SPDX-FileCopyrightText: 2026 Robert Marskar +// +// SPDX-License-Identifier: GPL-3.0-or-later + +// Benchmark: EBGeometry TriMeshSDF vs fcpw (https://github.com/rohan-sawhney/fcpw) on closest-point +// queries over a triangle mesh -- the task both libraries are built for (a BVH over triangles, +// answering "closest point on the surface"). Float precision (fcpw is float). +// +// The mesh is parsed once (shared, untimed) from the common-3d-test-models submodule. Each library +// then builds its own BVH over the same triangles (timed) and answers the same random closest-point +// queries (timed). Results are cross-checked: |TriMeshSDF signed distance| vs fcpw's unsigned +// closest-surface distance. Note TriMeshSDF additionally computes the sign (it is a signed-distance +// library); fcpw's findClosestPoint returns the unsigned closest point. + +#include +#include +#include +#include +#include +#include + +#include + +#include + +using T = float; // fcpw is float +using Vec3 = EBGeometry::Vec3T; +using Meta = EBGeometry::DCEL::DefaultMetaData; + +// SIMD-optimal branching factor / SoA triangle width for T (what readIntoTriangleBVH would pick). +constexpr std::size_t K = EBGeometry::BVH::DefaultBranchingRatio(); +constexpr std::size_t W = EBGeometry::TriangleSoA::DefaultWidth(); +using SDF = EBGeometry::TriMeshSDF; + +int +main(int argc, char** argv) +{ + const std::string objFile = (argc > 1) ? argv[1] : "../../Submodules/common-3d-test-models/data/armadillo.obj"; + constexpr std::size_t numQueries = 100000; + constexpr std::size_t sampleSize = 500; + constexpr std::size_t maxLeafGroups = 4; + + // Parse the mesh once -- shared preamble, not part of either library's timed build. + const auto tris = EBGeometry::Parser::readIntoTriangles(objFile); + const std::size_t nTri = tris.size(); + + std::printf("MeshSDF closest-point: EBGeometry TriMeshSDF vs fcpw\n"); + std::printf(" Mesh = %s (%zu triangles, float)\n Queries = %zu\n\n", objFile.c_str(), nTri, numQueries); + + EBGeometry::SimpleTimer timer; + + // Random query points in the mesh bounding box, expanded to 1.5x so queries sit inside and around. + Vec3 lo = +Vec3::max(); + Vec3 hi = -Vec3::max(); + for (const auto& tri : tris) { + for (const auto& p : tri->getVertexPositions()) { + lo = min(lo, p); + hi = max(hi, p); + } + } + const Vec3 center = T(0.5) * (lo + hi); + const Vec3 half = T(0.75) * (hi - lo); + std::mt19937 rng(12345u); + std::uniform_real_distribution u(T(-1), T(1)); + std::vector queries(numQueries); + for (auto& q : queries) { + q = center + Vec3(u(rng) * half[0], u(rng) * half[1], u(rng) * half[2]); + } + + // ── EBGeometry TriMeshSDF ── + timer.start(); + const SDF sdf(tris, EBGeometry::BVH::Build::SAH, maxLeafGroups); + timer.stop(); + const double ebBuildMs = 1.0e3 * timer.seconds(); + + std::vector ebDist(numQueries); + timer.start(); + for (std::size_t i = 0; i < numQueries; i++) { + ebDist[i] = std::abs(sdf.signedDistance(queries[i])); + } + timer.stop(); + const double ebQueryUs = 1.0e6 * timer.seconds() / double(numQueries); + + // ── fcpw ── + // fcpw needs vertices + triangle indices. Use an unwelded soup (3 vertices per triangle); closest + // point on triangles is unaffected by vertex sharing. Building this soup is fcpw's setup, not parse. + std::vector> V; + std::vector F; + V.reserve(3 * nTri); + F.reserve(nTri); + timer.start(); + for (const auto& tri : tris) { + const auto& p = tri->getVertexPositions(); + const int base = static_cast(V.size()); + for (int k = 0; k < 3; k++) { + V.emplace_back(fcpw::Vector<3>(p[k][0], p[k][1], p[k][2])); + } + F.emplace_back(fcpw::Vector3i(base, base + 1, base + 2)); + } + fcpw::Scene<3> scene; + scene.setObjectCount(1); + scene.setObjectVertices(V, 0); + scene.setObjectTriangles(F, 0); + scene.build(fcpw::AggregateType::Bvh_SurfaceArea, false /* vectorize (needs Enoki) */); + timer.stop(); + const double fcpwBuildMs = 1.0e3 * timer.seconds(); + + std::vector fcpwDist(numQueries); + timer.start(); + for (std::size_t i = 0; i < numQueries; i++) { + fcpw::Interaction<3> it; + scene.findClosestPoint(fcpw::Vector<3>(queries[i][0], queries[i][1], queries[i][2]), it); + fcpwDist[i] = it.d; + } + timer.stop(); + const double fcpwQueryUs = 1.0e6 * timer.seconds() / double(numQueries); + + // Cross-check on a spread sample: unsigned closest-surface distance must agree. + std::size_t bad = 0; + for (std::size_t s = 0; s < sampleSize; s++) { + const std::size_t i = s * (numQueries / sampleSize); + if (std::abs(ebDist[i] - fcpwDist[i]) > T(1.0e-3) * std::max(fcpwDist[i], T(1))) { + bad++; + } + } + + std::printf(" %-12s build %7.1f ms query %7.3f us/query\n", "TriMeshSDF", ebBuildMs, ebQueryUs); + std::printf(" %-12s build %7.1f ms query %7.3f us/query\n", "fcpw", fcpwBuildMs, fcpwQueryUs); + std::printf(" cross-check: %zu/%zu sample mismatches (|TriMeshSDF| vs fcpw unsigned distance)\n", bad, sampleSize); + return 0; +} diff --git a/Benchmark/NearestNeighbor/GNUmakefile b/Benchmark/NearestNeighbor/GNUmakefile new file mode 100644 index 00000000..b1236284 --- /dev/null +++ b/Benchmark/NearestNeighbor/GNUmakefile @@ -0,0 +1,24 @@ +# Benchmark: PointCloudBVH vs picoflann vs nanoflann (all-nearest-neighbor). +# Requires the Submodules/{picoflann,nanoflann} submodules: +# git submodule update --init Submodules/picoflann Submodules/nanoflann +# EBGEOMETRY_HOME defaults to the repo root (two levels up). + +EBGEOMETRY_HOME ?= ../.. +CXX ?= g++ +CXXFLAGS ?= -std=c++17 -O3 -march=native + +INCLUDES := -I$(EBGEOMETRY_HOME) \ + -I$(EBGEOMETRY_HOME)/Submodules/picoflann \ + -I$(EBGEOMETRY_HOME)/Submodules/nanoflann/include + +TARGET := NearestNeighbor.ex + +$(TARGET): main.cpp + $(CXX) $(CXXFLAGS) $(INCLUDES) $< -o $@ + +.PHONY: run clean +run: $(TARGET) + ./$(TARGET) + +clean: + $(RM) $(TARGET) diff --git a/Benchmark/NearestNeighbor/main.cpp b/Benchmark/NearestNeighbor/main.cpp new file mode 100644 index 00000000..9b1bcb50 --- /dev/null +++ b/Benchmark/NearestNeighbor/main.cpp @@ -0,0 +1,222 @@ +// SPDX-FileCopyrightText: 2026 Robert Marskar +// +// SPDX-License-Identifier: GPL-3.0-or-later + +// Benchmark: EBGeometry PointCloudBVH vs picoflann vs nanoflann, all-nearest-neighbor. +// 500k random 3D points in the unit cube (double). Every point's nearest OTHER point. All three +// verified against a brute-force sample. Both KD-trees used vanilla; distances are squared throughout. + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "nanoflann.hpp" +#include "picoflann.h" + +using T = double; +using Vec3 = EBGeometry::Vec3T; + +constexpr std::size_t numPoints = 500000; +constexpr std::size_t sampleSize = 500; +constexpr std::uint64_t pointSeed = 123456789ULL; + +namespace { + +// picoflann adapter. +struct Vec3Adapter +{ + inline T + operator()(const Vec3& a_p, int a_dim) const + { + return a_p[a_dim]; + } +}; + +// nanoflann dataset adaptor over std::vector. +struct NanoCloud +{ + const std::vector& pts; + inline std::size_t + kdtree_get_point_count() const + { + return pts.size(); + } + inline T + kdtree_get_pt(const std::size_t a_idx, const std::size_t a_dim) const + { + return pts[a_idx][a_dim]; + } + template + bool + kdtree_get_bbox(BBOX&) const + { + return false; + } +}; +using NanoTree = nanoflann::KDTreeSingleIndexAdaptor, NanoCloud, 3>; + +T +bruteForceNN2(std::size_t a_self, const std::vector& a_pos) +{ + T best = std::numeric_limits::max(); + for (std::size_t i = 0; i < a_pos.size(); i++) { + if (i != a_self) { + best = std::min(best, (a_pos[i] - a_pos[a_self]).length2()); + } + } + return best; +} + +} // namespace + +int +main() +{ + std::cout << "All-nearest-neighbor: PointCloudBVH vs picoflann vs nanoflann\n"; + std::cout << " Points = " << numPoints << " (double, unit cube)\n\n"; + + const std::vector positions = EBGeometry::Random::samplePoints(numPoints, pointSeed); + const std::vector meta(numPoints); + + EBGeometry::SimpleTimer timer; + + // Query the flann trees in a spatially-coherent (Hilbert) order too, so their node cache is as warm + // as EBGeometry's leaf-order batch. EBGeometry gets its order free from the build; the flann libs + // must sort -- time that once so it can be folded in if desired. + timer.start(); + const std::vector order = EBGeometry::SFC::order(positions); + timer.stop(); + const double sortUsPerPt = 1.0e6 * timer.seconds() / double(numPoints); + + const std::size_t stride = numPoints / sampleSize; + std::vector truth(sampleSize); + timer.start(); + for (std::size_t s = 0; s < sampleSize; s++) { + truth[s] = bruteForceNN2(s * stride, positions); + } + timer.stop(); + const double bruteUsPerPt = 1.0e6 * timer.seconds() / double(sampleSize); + + auto ok = [&](T a_got, std::size_t a_s) { + return std::abs(a_got - truth[a_s]) <= 1.0e-9 * std::max(truth[a_s], T(1)); + }; + + std::cout << std::left << std::setw(22) << "Method" << std::right << std::setw(12) << "Build(ms)" << std::setw(14) + << "Query(us/pt)" << std::setw(12) << "vs brute" << '\n'; + std::cout << std::string(60, '-') << '\n'; + std::cout << std::fixed; + std::cout << std::left << std::setw(22) << "Brute force" << std::right << std::setw(12) << "--" << std::setw(14) + << std::setprecision(3) << bruteUsPerPt << std::setw(12) << "1.0x" << '\n'; + + auto row = [&](const char* a_name, double a_buildMs, double a_queryUs, std::size_t a_bad) { + std::cout << std::left << std::setw(22) << a_name << std::right << std::setw(12) << std::setprecision(1) + << a_buildMs << std::setw(14) << std::setprecision(3) << a_queryUs << std::setw(11) + << std::setprecision(1) << bruteUsPerPt / a_queryUs << "x" << (a_bad ? " MISMATCH!" : "") << '\n'; + }; + + // ── EBGeometry PointCloudBVH (batched all-NN) ── + { + timer.start(); + const EBGeometry::PointCloudBVH bvh(positions, meta); + timer.stop(); + const double buildMs = 1.0e3 * timer.seconds(); + + timer.start(); + const auto graph = bvh.allNearestNeighbors(1); + timer.stop(); + const double queryUs = 1.0e6 * timer.seconds() / double(numPoints); + + std::size_t bad = 0; + for (std::size_t s = 0; s < sampleSize; s++) { + bad += !ok(graph[s * stride].distanceSquared, s); + } + row("PointCloudBVH", buildMs, queryUs, bad); + } + + // ── picoflann (per-point searchKnn) ── + { + picoflann::KdTreeIndex<3, Vec3Adapter> kdtree; + timer.start(); + kdtree.build(positions); + timer.stop(); + const double buildMs = 1.0e3 * timer.seconds(); + + volatile T sink = T(0); + timer.start(); + for (const std::uint32_t p : order) { + const auto res = kdtree.searchKnn(positions, positions[p], 2); + for (const auto& pr : res) { + if (pr.first != p) { + sink += pr.second; + break; + } + } + } + timer.stop(); + (void)sink; + const double queryUs = 1.0e6 * timer.seconds() / double(numPoints); + + std::size_t bad = 0; + for (std::size_t s = 0; s < sampleSize; s++) { + const auto res = kdtree.searchKnn(positions, positions[s * stride], 2); + T got = std::numeric_limits::max(); + for (const auto& pr : res) { + if (pr.first != s * stride) { + got = pr.second; + break; + } + } + bad += !ok(got, s); + } + row("picoflann", buildMs, queryUs, bad); + } + + // ── nanoflann (per-point findNeighbors, k=2) ── + { + NanoCloud cloud{positions}; + timer.start(); + NanoTree index(3, cloud, nanoflann::KDTreeSingleIndexAdaptorParams(10 /* leaf_max_size */)); + index.buildIndex(); + timer.stop(); + const double buildMs = 1.0e3 * timer.seconds(); + + auto nnOther = [&](std::size_t i) { + std::size_t idx[2]; + T d2[2]; + nanoflann::KNNResultSet rs(2); + rs.init(idx, d2); + const T qp[3] = {positions[i][0], positions[i][1], positions[i][2]}; + index.findNeighbors(rs, qp); + // idx[0] is the point itself (distance 0); take the first neighbor that is not itself. + return (idx[0] != i) ? d2[0] : d2[1]; + }; + + volatile T sink = T(0); + timer.start(); + for (const std::uint32_t p : order) { + sink += nnOther(p); + } + timer.stop(); + (void)sink; + const double queryUs = 1.0e6 * timer.seconds() / double(numPoints); + + std::size_t bad = 0; + for (std::size_t s = 0; s < sampleSize; s++) { + bad += !ok(nnOther(s * stride), s); + } + row("nanoflann", buildMs, queryUs, bad); + } + + std::cout << "\n flann query loops iterate in Hilbert order (warm cache, like EBGeometry's leaf order).\n"; + std::cout << " One-time Hilbert sort the flann libs need for that order: " << std::setprecision(3) << sortUsPerPt + << " us/pt (add to their query if counted;\n EBGeometry reuses its build order for free).\n"; + + return 0; +} diff --git a/Benchmark/README.md b/Benchmark/README.md new file mode 100644 index 00000000..2c04bf86 --- /dev/null +++ b/Benchmark/README.md @@ -0,0 +1,67 @@ +Benchmark +========= + +Head-to-head performance comparisons of EBGeometry against other open-source geometry-query +libraries. Like `Integrations/`, this is **illustrative and not part of CI** — it depends on +external libraries with their own build systems and release cadences (see issue #109). + +The comparison libraries are pinned as git submodules under [`Submodules/`](../Submodules): + +* [nanoflann](https://github.com/jlblancoc/nanoflann) — header-only KD-tree (point kNN) +* [picoflann](https://github.com/rmsalinas/picoflann) — tiny header-only KD-tree (point kNN) +* [fcpw](https://github.com/rohan-sawhney/fcpw) — closest-point / SDF on triangle meshes + +Fetch them (and the mesh submodule) with: + +```bash +git submodule update --init Submodules/nanoflann Submodules/picoflann Submodules/fcpw Submodules/common-3d-test-models +git -C Submodules/fcpw submodule update --init deps/eigen # fcpw's Eigen (skip the GPU slang-rhi dep) +``` + +Each benchmark has a `GNUmakefile` (`make && ./.ex`). Every result is cross-checked against a +brute-force / independent baseline so a wrong answer shows up as a mismatch. + +`NearestNeighbor/` — all-nearest-neighbor on a point cloud +---------------------------------------------------------- + +For every point in a 500,000-point uniform cloud in the unit cube (double precision), find its +nearest *other* point. Compares `PointCloudBVH` vs picoflann vs nanoflann. + +Representative result (one machine): + +``` +Method Build(ms) Query(us/pt) +PointCloudBVH ~78 0.29 +nanoflann ~120 0.25 +picoflann ~65 0.47 +``` + +- The KD-tree queries are iterated in **Hilbert order** so their node cache is as warm as + `PointCloudBVH`'s leaf-order batch (natural order is ~2x slower — unfair). That order costs the + KD-trees a one-time ~0.18 us/pt spatial sort; `PointCloudBVH` reuses the ordering its build already + produced, for free. Folding the sort in: `PointCloudBVH` ~0.29 vs nanoflann ~0.43 vs picoflann ~0.65. +- nanoflann has the fastest raw traversal; `PointCloudBVH` wins the end-to-end all-NN job (free query + order + faster build) and builds ~1.5x faster than nanoflann. + +`MeshSDF/` — closest-point on a triangle mesh +--------------------------------------------- + +Closest-point queries against a triangle mesh (armadillo, ~100k triangles, float) — the task fcpw and +EBGeometry's `TriMeshSDF` are both built for. The mesh is parsed once; each library then builds its +own BVH over the same triangles and answers the same queries. + +Representative result (one machine, 100k queries): + +``` +Method Build(ms) Query(us/query) +TriMeshSDF ~65 2.6 +fcpw ~55 7.5 +``` + +- **Caveat:** fcpw here runs its **scalar Eigen fallback** (`FCPW_USE_ENOKI` off), while `TriMeshSDF` + is SIMD-vectorized (SoA triangle leaves). So the query gap flatters EBGeometry — fcpw's intended + fast path uses Enoki CPU vectorization, which would narrow it. Enabling Enoki is a heavier build + (an extra dependency) and is left as a follow-up. Build times and correctness are directly + comparable; the cross-check confirms the distances agree. +- `TriMeshSDF` computes the *signed* distance (its purpose); fcpw's `findClosestPoint` returns the + unsigned closest point, so the comparison is on unsigned closest-surface distance. diff --git a/Docs/Sphinx/source/Benchmark.rst b/Docs/Sphinx/source/Benchmark.rst new file mode 100644 index 00000000..5115a975 --- /dev/null +++ b/Docs/Sphinx/source/Benchmark.rst @@ -0,0 +1,29 @@ +.. _Chap:Benchmark: + +Benchmarks +========== + +.. important:: + + Like the :ref:`Chap:Integrations` examples, the benchmarks are illustrative and **not built or run + as part of EBGeometry's continuous integration** -- they depend on external libraries with their + own build systems. They pin those libraries as git submodules under :file:`Submodules/`. See the + benchmark tracking issue on GitHub for context and planned additions. + +The :file:`Benchmark/` folder compares EBGeometry against other open-source geometry-query libraries +on tasks they have in common. Fetch the comparison libraries (and the mesh submodule) with: + +.. code-block:: bash + + git submodule update --init Submodules/nanoflann Submodules/picoflann Submodules/fcpw Submodules/common-3d-test-models + git -C Submodules/fcpw submodule update --init deps/eigen # fcpw's Eigen (skip the GPU dep) + +Each benchmark ships a ``GNUmakefile`` (``make && ./.ex``) and cross-checks every result against +a brute-force / independent baseline. See each folder's ``README.md`` for the full detail and +representative numbers. + +* :file:`Benchmark/NearestNeighbor` -- all-nearest-neighbor over a point cloud: ``PointCloudBVH`` vs + `nanoflann `_ vs + `picoflann `_. +* :file:`Benchmark/MeshSDF` -- closest-point on a triangle mesh: ``TriMeshSDF`` vs + `fcpw `_. diff --git a/Docs/Sphinx/source/index.rst b/Docs/Sphinx/source/index.rst index b54838ee..d4d49063 100644 --- a/Docs/Sphinx/source/index.rst +++ b/Docs/Sphinx/source/index.rst @@ -128,6 +128,7 @@ Examples ExampleNearestNeighborBVH.rst ExampleNearestNeighborHashGrid.rst Integrations.rst + Benchmark.rst Contributing and testing ************************ diff --git a/REUSE.toml b/REUSE.toml index 9f53e573..08c965ab 100644 --- a/REUSE.toml +++ b/REUSE.toml @@ -69,6 +69,16 @@ path = [ SPDX-FileCopyrightText = "2022 Robert Marskar " SPDX-License-Identifier = "GPL-3.0-or-later" +# Benchmark comparison programs (illustrative, not CI-tested). main.cpp files carry inline SPDX. +[[annotations]] +path = [ + "Benchmark/README.md", + "Benchmark/**/GNUmakefile", + "Benchmark/**/README.md", +] +SPDX-FileCopyrightText = "2022 Robert Marskar " +SPDX-License-Identifier = "GPL-3.0-or-later" + # Tests: build files and the small self-authored test fixture mesh [[annotations]] path = [ From 1d7bf72040a69050af9d685d9f56118d7a5f3a77 Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 12 Jul 2026 01:36:15 +0200 Subject: [PATCH 03/17] Benchmark/MeshSDF: build fcpw with Enoki and add TriangleMeshDistance to the comparison - fcpw is now built with its Enoki CPU vectorization (FCPW_USE_ENOKI, vectorized MBVH) instead of the scalar Eigen fallback -- a like-for-like float/SIMD comparison with TriMeshSDF. Enoki's headers are vendored inside the fcpw submodule, so no extra fetch is needed; FCPW_SIMD_WIDTH is a GNUmakefile variable (default 8 = AVX2; set 4 for SSE, 16 for AVX-512). - Add TriangleMeshDistance (https://github.com/InteractiveComputerGraphics/TriangleMeshDistance) as a third contender, pinned as a submodule under Submodules/. It is a header-only, double-precision, scalar signed-distance library -- a widely-used reference point, noted as not a same-precision comparison. Its distances are cross-checked against TriMeshSDF alongside fcpw's (0 mismatches). Representative (armadillo ~100k triangles, 100k queries, one machine): TriMeshSDF build ~63 ms query ~2.6 us fcpw (Enoki) build ~46 ms query ~3.5 us (was ~7.5 us scalar) TriangleMeshDistance build ~115 ms query ~10.5 us (double, scalar) Updates the Benchmark README and Benchmark.rst accordingly. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01HS5LFQihKoMrE6nbPaQW9R --- .gitmodules | 3 ++ Benchmark/MeshSDF/GNUmakefile | 19 +++++--- Benchmark/MeshSDF/main.cpp | 79 ++++++++++++++++++++++++++------ Benchmark/README.md | 36 +++++++++------ Docs/Sphinx/source/Benchmark.rst | 6 ++- Submodules/TriangleMeshDistance | 1 + 6 files changed, 108 insertions(+), 36 deletions(-) create mode 160000 Submodules/TriangleMeshDistance diff --git a/.gitmodules b/.gitmodules index ed5eb559..061ead1a 100644 --- a/.gitmodules +++ b/.gitmodules @@ -14,3 +14,6 @@ path = Submodules/fcpw url = https://github.com/rohan-sawhney/fcpw.git shallow = true +[submodule "Submodules/TriangleMeshDistance"] + path = Submodules/TriangleMeshDistance + url = https://github.com/InteractiveComputerGraphics/TriangleMeshDistance diff --git a/Benchmark/MeshSDF/GNUmakefile b/Benchmark/MeshSDF/GNUmakefile index 8ec70447..712f1768 100644 --- a/Benchmark/MeshSDF/GNUmakefile +++ b/Benchmark/MeshSDF/GNUmakefile @@ -1,21 +1,28 @@ -# Benchmark: EBGeometry TriMeshSDF vs fcpw (closest-point on a triangle mesh). -# Requires the Submodules/fcpw submodule plus its Eigen dependency, and the mesh submodule: -# git submodule update --init Submodules/fcpw Submodules/common-3d-test-models +# Benchmark: EBGeometry TriMeshSDF vs fcpw vs TriangleMeshDistance (closest-point on a triangle mesh). +# Requires the mesh submodule, fcpw (Enoki is vendored inside fcpw; Eigen is a nested submodule), and +# TriangleMeshDistance: +# git submodule update --init Submodules/fcpw Submodules/common-3d-test-models Submodules/TriangleMeshDistance # git -C Submodules/fcpw submodule update --init deps/eigen # Eigen only (skip GPU slang-rhi) -# EBGEOMETRY_HOME defaults to the repo root (two levels up). +# fcpw is built with its Enoki CPU vectorization: FCPW_SIMD_WIDTH matches the ISA (4=SSE, 8=AVX2, +# 16=AVX-512); adjust it below to your machine. EBGEOMETRY_HOME defaults to the repo root. EBGEOMETRY_HOME ?= ../.. CXX ?= g++ CXXFLAGS ?= -std=c++17 -O3 -march=native +FCPW_SIMD_WIDTH ?= 8 INCLUDES := -I$(EBGEOMETRY_HOME) \ -I$(EBGEOMETRY_HOME)/Submodules/fcpw/include \ - -I$(EBGEOMETRY_HOME)/Submodules/fcpw/deps/eigen + -I$(EBGEOMETRY_HOME)/Submodules/fcpw/deps/eigen \ + -I$(EBGEOMETRY_HOME)/Submodules/fcpw/deps/enoki/include \ + -I$(EBGEOMETRY_HOME)/Submodules/TriangleMeshDistance/TriangleMeshDistance/include + +DEFINES := -DFCPW_USE_ENOKI -DFCPW_SIMD_WIDTH=$(FCPW_SIMD_WIDTH) TARGET := MeshSDF.ex $(TARGET): main.cpp - $(CXX) $(CXXFLAGS) $(INCLUDES) $< -o $@ + $(CXX) $(CXXFLAGS) $(DEFINES) $(INCLUDES) $< -o $@ .PHONY: run clean run: $(TARGET) diff --git a/Benchmark/MeshSDF/main.cpp b/Benchmark/MeshSDF/main.cpp index 02713739..b3dd0c3d 100644 --- a/Benchmark/MeshSDF/main.cpp +++ b/Benchmark/MeshSDF/main.cpp @@ -2,15 +2,21 @@ // // SPDX-License-Identifier: GPL-3.0-or-later -// Benchmark: EBGeometry TriMeshSDF vs fcpw (https://github.com/rohan-sawhney/fcpw) on closest-point -// queries over a triangle mesh -- the task both libraries are built for (a BVH over triangles, -// answering "closest point on the surface"). Float precision (fcpw is float). +// Benchmark: EBGeometry TriMeshSDF vs fcpw (https://github.com/rohan-sawhney/fcpw) vs +// TriangleMeshDistance (https://github.com/InteractiveComputerGraphics/TriangleMeshDistance) on +// closest-point queries over a triangle mesh -- the task all three are built for (an acceleration +// structure over triangles, answering "closest point on the surface"). // // The mesh is parsed once (shared, untimed) from the common-3d-test-models submodule. Each library -// then builds its own BVH over the same triangles (timed) and answers the same random closest-point -// queries (timed). Results are cross-checked: |TriMeshSDF signed distance| vs fcpw's unsigned -// closest-surface distance. Note TriMeshSDF additionally computes the sign (it is a signed-distance -// library); fcpw's findClosestPoint returns the unsigned closest point. +// then builds its own structure over the same triangles (timed) and answers the same random +// closest-point queries (timed). Results are cross-checked against TriMeshSDF's unsigned distance. +// +// Precision / vectorization caveats (each library on its own intended fast path): +// - TriMeshSDF -- float, SIMD-vectorized (this build's native ISA). +// - fcpw -- float; built with its Enoki CPU vectorization (FCPW_USE_ENOKI, vectorized BVH). +// fcpw returns an unsigned closest point; TriMeshSDF additionally computes a sign. +// - TriangleMeshDistance -- double, scalar (header-only, no SIMD). A signed-distance library like +// TriMeshSDF; runs its queries in double, so it is not a same-precision comparison. #include #include @@ -22,6 +28,7 @@ #include #include +#include using T = float; // fcpw is float using Vec3 = EBGeometry::Vec3T; @@ -101,7 +108,7 @@ main(int argc, char** argv) scene.setObjectCount(1); scene.setObjectVertices(V, 0); scene.setObjectTriangles(F, 0); - scene.build(fcpw::AggregateType::Bvh_SurfaceArea, false /* vectorize (needs Enoki) */); + scene.build(fcpw::AggregateType::Bvh_SurfaceArea, true /* vectorize (Enoki MBVH) */); timer.stop(); const double fcpwBuildMs = 1.0e3 * timer.seconds(); @@ -115,17 +122,61 @@ main(int argc, char** argv) timer.stop(); const double fcpwQueryUs = 1.0e6 * timer.seconds() / double(numQueries); - // Cross-check on a spread sample: unsigned closest-surface distance must agree. - std::size_t bad = 0; + // ── TriangleMeshDistance (double, scalar, header-only) ── + // Same unwelded triangle soup, in double. Building the arrays + the structure is its timed setup, + // mirroring how fcpw's soup construction is folded into fcpw's build above. + std::vector tmdVertices; + std::vector tmdTriangles; + tmdVertices.reserve(9 * nTri); + tmdTriangles.reserve(3 * nTri); + timer.start(); + for (const auto& tri : tris) { + const auto& p = tri->getVertexPositions(); + const int base = static_cast(tmdVertices.size() / 3); + for (int k = 0; k < 3; k++) { + tmdVertices.push_back(double(p[k][0])); + tmdVertices.push_back(double(p[k][1])); + tmdVertices.push_back(double(p[k][2])); + } + tmdTriangles.push_back(base); + tmdTriangles.push_back(base + 1); + tmdTriangles.push_back(base + 2); + } + const tmd::TriangleMeshDistance tmdMesh( + tmdVertices.data(), tmdVertices.size() / 3, tmdTriangles.data(), tmdTriangles.size() / 3); + timer.stop(); + const double tmdBuildMs = 1.0e3 * timer.seconds(); + + std::vector tmdDist(numQueries); + timer.start(); + for (std::size_t i = 0; i < numQueries; i++) { + const tmd::Result r = + tmdMesh.signed_distance({double(queries[i][0]), double(queries[i][1]), double(queries[i][2])}); + tmdDist[i] = T(std::abs(r.distance)); + } + timer.stop(); + const double tmdQueryUs = 1.0e6 * timer.seconds() / double(numQueries); + + // Cross-check on a spread sample: unsigned closest-surface distance must agree across all three. + std::size_t badFcpw = 0; + std::size_t badTmd = 0; for (std::size_t s = 0; s < sampleSize; s++) { const std::size_t i = s * (numQueries / sampleSize); if (std::abs(ebDist[i] - fcpwDist[i]) > T(1.0e-3) * std::max(fcpwDist[i], T(1))) { - bad++; + badFcpw++; + } + if (std::abs(ebDist[i] - tmdDist[i]) > T(1.0e-3) * std::max(tmdDist[i], T(1))) { + badTmd++; } } - std::printf(" %-12s build %7.1f ms query %7.3f us/query\n", "TriMeshSDF", ebBuildMs, ebQueryUs); - std::printf(" %-12s build %7.1f ms query %7.3f us/query\n", "fcpw", fcpwBuildMs, fcpwQueryUs); - std::printf(" cross-check: %zu/%zu sample mismatches (|TriMeshSDF| vs fcpw unsigned distance)\n", bad, sampleSize); + std::printf(" %-22s build %7.1f ms query %7.3f us/query\n", "TriMeshSDF", ebBuildMs, ebQueryUs); + std::printf(" %-22s build %7.1f ms query %7.3f us/query\n", "fcpw (Enoki)", fcpwBuildMs, fcpwQueryUs); + std::printf(" %-22s build %7.1f ms query %7.3f us/query\n", "TriangleMeshDistance", tmdBuildMs, tmdQueryUs); + std::printf(" cross-check vs TriMeshSDF: fcpw %zu/%zu, TriangleMeshDistance %zu/%zu sample mismatches\n", + badFcpw, + sampleSize, + badTmd, + sampleSize); return 0; } diff --git a/Benchmark/README.md b/Benchmark/README.md index 2c04bf86..47aa5e0d 100644 --- a/Benchmark/README.md +++ b/Benchmark/README.md @@ -10,14 +10,19 @@ The comparison libraries are pinned as git submodules under [`Submodules/`](../S * [nanoflann](https://github.com/jlblancoc/nanoflann) — header-only KD-tree (point kNN) * [picoflann](https://github.com/rmsalinas/picoflann) — tiny header-only KD-tree (point kNN) * [fcpw](https://github.com/rohan-sawhney/fcpw) — closest-point / SDF on triangle meshes +* [TriangleMeshDistance](https://github.com/InteractiveComputerGraphics/TriangleMeshDistance) — header-only signed distance to triangle meshes Fetch them (and the mesh submodule) with: ```bash -git submodule update --init Submodules/nanoflann Submodules/picoflann Submodules/fcpw Submodules/common-3d-test-models +git submodule update --init Submodules/nanoflann Submodules/picoflann Submodules/fcpw \ + Submodules/TriangleMeshDistance Submodules/common-3d-test-models git -C Submodules/fcpw submodule update --init deps/eigen # fcpw's Eigen (skip the GPU slang-rhi dep) ``` +fcpw's Enoki CPU-vectorization headers are vendored inside the fcpw submodule, so no extra fetch is +needed for it. + Each benchmark has a `GNUmakefile` (`make && ./.ex`). Every result is cross-checked against a brute-force / independent baseline so a wrong answer shows up as a mismatch. @@ -46,22 +51,25 @@ picoflann ~65 0.47 `MeshSDF/` — closest-point on a triangle mesh --------------------------------------------- -Closest-point queries against a triangle mesh (armadillo, ~100k triangles, float) — the task fcpw and -EBGeometry's `TriMeshSDF` are both built for. The mesh is parsed once; each library then builds its -own BVH over the same triangles and answers the same queries. +Closest-point queries against a triangle mesh (armadillo, ~100k triangles) — the task `TriMeshSDF`, +fcpw, and TriangleMeshDistance are all built for. The mesh is parsed once; each library then builds +its own structure over the same triangles and answers the same queries. Representative result (one machine, 100k queries): ``` -Method Build(ms) Query(us/query) -TriMeshSDF ~65 2.6 -fcpw ~55 7.5 +Method Build(ms) Query(us/query) +TriMeshSDF ~63 2.6 +fcpw (Enoki) ~46 3.5 +TriangleMeshDistance ~115 10.5 ``` -- **Caveat:** fcpw here runs its **scalar Eigen fallback** (`FCPW_USE_ENOKI` off), while `TriMeshSDF` - is SIMD-vectorized (SoA triangle leaves). So the query gap flatters EBGeometry — fcpw's intended - fast path uses Enoki CPU vectorization, which would narrow it. Enabling Enoki is a heavier build - (an extra dependency) and is left as a follow-up. Build times and correctness are directly - comparable; the cross-check confirms the distances agree. -- `TriMeshSDF` computes the *signed* distance (its purpose); fcpw's `findClosestPoint` returns the - unsigned closest point, so the comparison is on unsigned closest-surface distance. +- **fcpw** is built with its **Enoki CPU vectorization** (`FCPW_USE_ENOKI`, vectorized MBVH), the same + SIMD fast path `TriMeshSDF` uses — a like-for-like float/SIMD comparison. (`FCPW_SIMD_WIDTH` in the + `GNUmakefile` should match your ISA: 4=SSE, 8=AVX2, 16=AVX-512.) fcpw's `findClosestPoint` returns + the *unsigned* closest point. +- **TriangleMeshDistance** is a header-only, **double-precision, scalar** (non-SIMD) signed-distance + library — so its query runs in double and is not a same-precision comparison; it is included as a + widely-used point of reference. +- `TriMeshSDF` and TriangleMeshDistance compute the *signed* distance (their purpose); the comparison + is on unsigned closest-surface distance, and the cross-check confirms all three agree. diff --git a/Docs/Sphinx/source/Benchmark.rst b/Docs/Sphinx/source/Benchmark.rst index 5115a975..fed6ae7d 100644 --- a/Docs/Sphinx/source/Benchmark.rst +++ b/Docs/Sphinx/source/Benchmark.rst @@ -15,7 +15,8 @@ on tasks they have in common. Fetch the comparison libraries (and the mesh submo .. code-block:: bash - git submodule update --init Submodules/nanoflann Submodules/picoflann Submodules/fcpw Submodules/common-3d-test-models + git submodule update --init Submodules/nanoflann Submodules/picoflann Submodules/fcpw \ + Submodules/TriangleMeshDistance Submodules/common-3d-test-models git -C Submodules/fcpw submodule update --init deps/eigen # fcpw's Eigen (skip the GPU dep) Each benchmark ships a ``GNUmakefile`` (``make && ./.ex``) and cross-checks every result against @@ -26,4 +27,5 @@ representative numbers. `nanoflann `_ vs `picoflann `_. * :file:`Benchmark/MeshSDF` -- closest-point on a triangle mesh: ``TriMeshSDF`` vs - `fcpw `_. + `fcpw `_ (built with its Enoki CPU vectorization) vs + `TriangleMeshDistance `_. diff --git a/Submodules/TriangleMeshDistance b/Submodules/TriangleMeshDistance new file mode 160000 index 00000000..5530eaa8 --- /dev/null +++ b/Submodules/TriangleMeshDistance @@ -0,0 +1 @@ +Subproject commit 5530eaa85537cc62b884c3f8a3e64e000129eb93 From e1c355de55d11a2a9eca5ba2cc86798eed952e1b Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 15 Jul 2026 19:55:19 +0000 Subject: [PATCH 04/17] Relocate benchmark submodules under Benchmark/, soften framing, add compile-only CI Move the comparison-library submodules (nanoflann, picoflann, fcpw, TriangleMeshDistance) out of Submodules/ to live directly under Benchmark/, and move common-3d-test-models back to the repository top level where it has always been. Update .gitmodules, REUSE.toml, all example/integration path references, and the benchmark docs to match. Soften the benchmark framing in Benchmark/README.md and Docs/Sphinx/source/Benchmark.rst: present the suite as reproducible, cross-checked methodology whose single-machine numbers are illustrative snapshots rather than a scoreboard. Add a compile-only Benchmark-Compile CI job that fetches each benchmark's comparison-library submodules and builds (but never runs or times) the two programs via their GNUmakefiles, as a smoke test against bit-rot. Document the new job and the softened framing in ContinuousIntegration.rst. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01JhdRJydutVJqQjg2Fuee9Z --- .github/workflows/CI.yml | 33 ++++++++++++++++- .gitmodules | 18 +++++----- Benchmark/MeshSDF/GNUmakefile | 12 +++---- Benchmark/MeshSDF/main.cpp | 2 +- Benchmark/NearestNeighbor/GNUmakefile | 8 ++--- Benchmark/README.md | 36 +++++++++++-------- .../TriangleMeshDistance | 0 {Submodules => Benchmark}/fcpw | 0 {Submodules => Benchmark}/nanoflann | 0 {Submodules => Benchmark}/picoflann | 0 Docs/Sphinx/source/Benchmark.rst | 30 +++++++++------- Docs/Sphinx/source/ContinuousIntegration.rst | 19 +++++++--- Docs/Sphinx/source/ExampleCSGUnion.rst | 2 +- Docs/Sphinx/source/ExampleMeshSDF.rst | 2 +- Docs/Sphinx/source/ObtainingEBGeometry.rst | 6 ++-- Examples/CMakeLists.txt | 2 +- Examples/CSGUnion/README.md | 2 +- Examples/CSGUnion/main.cpp | 2 +- Examples/MeshSDF/README.md | 2 +- Examples/MeshSDF/main.cpp | 4 +-- Examples/README.md | 2 +- Integrations/AMReX/MeshSDF/main.cpp | 2 +- Integrations/AMReX/PaintEB/main.cpp | 2 +- REUSE.toml | 2 +- ...on-3d-test-models => common-3d-test-models | 0 25 files changed, 121 insertions(+), 67 deletions(-) rename {Submodules => Benchmark}/TriangleMeshDistance (100%) rename {Submodules => Benchmark}/fcpw (100%) rename {Submodules => Benchmark}/nanoflann (100%) rename {Submodules => Benchmark}/picoflann (100%) rename Submodules/common-3d-test-models => common-3d-test-models (100%) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index e4291c2e..41c5c3f6 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -430,8 +430,39 @@ jobs: UBSAN_OPTIONS: print_stacktrace=1:halt_on_error=1 ASAN_OPTIONS: detect_leaks=1 + Benchmark-Compile: + needs: [Formatting, Codespell, Reuse, Doxygen-check] + name: Benchmarks (compile-only, ${{ matrix.directory }}) + runs-on: ubuntu-latest + timeout-minutes: 20 + # The benchmarks link against external comparison libraries with their own build + # systems, so they are not part of the normal test build. This job only *compiles* + # them (a smoke test against bit-rot as the library evolves) -- it never runs or + # times them, so no mesh data (common-3d-test-models) is needed and results are + # never treated as measurements. See Benchmark/README.md. + strategy: + fail-fast: false + matrix: + include: + - directory: Benchmark/NearestNeighbor + submodules: Benchmark/picoflann Benchmark/nanoflann + - directory: Benchmark/MeshSDF + submodules: Benchmark/fcpw Benchmark/TriangleMeshDistance + steps: + - uses: actions/checkout@v7 + - name: Fetch comparison-library submodules + run: | + git submodule update --init --depth 1 ${{ matrix.submodules }} + if [ "${{ matrix.directory }}" = "Benchmark/MeshSDF" ]; then + # fcpw's Eigen (nested submodule); Enoki is vendored inside fcpw, GPU slang-rhi is skipped. + git -C Benchmark/fcpw submodule update --init --depth 1 deps/eigen + fi + - name: Compile benchmark (build only, do not run) + working-directory: ${{ matrix.directory }} + run: make + CI-passed: - needs: [Formatting, Codespell, Reuse, Doxygen-check, Linux-GNU, Linux-Intel, Examples-GNUMake, Examples-CMake, Examples-FloatPrecision, Build-documentation, Unit-Tests, Release-Test, Sanitizers] + needs: [Formatting, Codespell, Reuse, Doxygen-check, Linux-GNU, Linux-Intel, Examples-GNUMake, Examples-CMake, Examples-FloatPrecision, Build-documentation, Unit-Tests, Release-Test, Sanitizers, Benchmark-Compile] runs-on: ubuntu-latest steps: - name: Do nothing diff --git a/.gitmodules b/.gitmodules index 061ead1a..8a654ef1 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,19 +1,19 @@ [submodule "common-3d-test-models"] - path = Submodules/common-3d-test-models + path = common-3d-test-models url = https://github.com/alecjacobson/common-3d-test-models shallow = true -[submodule "Submodules/nanoflann"] - path = Submodules/nanoflann +[submodule "Benchmark/nanoflann"] + path = Benchmark/nanoflann url = https://github.com/jlblancoc/nanoflann.git shallow = true -[submodule "Submodules/picoflann"] - path = Submodules/picoflann +[submodule "Benchmark/picoflann"] + path = Benchmark/picoflann url = https://github.com/rmsalinas/picoflann.git shallow = true -[submodule "Submodules/fcpw"] - path = Submodules/fcpw +[submodule "Benchmark/fcpw"] + path = Benchmark/fcpw url = https://github.com/rohan-sawhney/fcpw.git shallow = true -[submodule "Submodules/TriangleMeshDistance"] - path = Submodules/TriangleMeshDistance +[submodule "Benchmark/TriangleMeshDistance"] + path = Benchmark/TriangleMeshDistance url = https://github.com/InteractiveComputerGraphics/TriangleMeshDistance diff --git a/Benchmark/MeshSDF/GNUmakefile b/Benchmark/MeshSDF/GNUmakefile index 712f1768..b6150a1c 100644 --- a/Benchmark/MeshSDF/GNUmakefile +++ b/Benchmark/MeshSDF/GNUmakefile @@ -1,8 +1,8 @@ # Benchmark: EBGeometry TriMeshSDF vs fcpw vs TriangleMeshDistance (closest-point on a triangle mesh). # Requires the mesh submodule, fcpw (Enoki is vendored inside fcpw; Eigen is a nested submodule), and # TriangleMeshDistance: -# git submodule update --init Submodules/fcpw Submodules/common-3d-test-models Submodules/TriangleMeshDistance -# git -C Submodules/fcpw submodule update --init deps/eigen # Eigen only (skip GPU slang-rhi) +# git submodule update --init Benchmark/fcpw common-3d-test-models Benchmark/TriangleMeshDistance +# git -C Benchmark/fcpw submodule update --init deps/eigen # Eigen only (skip GPU slang-rhi) # fcpw is built with its Enoki CPU vectorization: FCPW_SIMD_WIDTH matches the ISA (4=SSE, 8=AVX2, # 16=AVX-512); adjust it below to your machine. EBGEOMETRY_HOME defaults to the repo root. @@ -12,10 +12,10 @@ CXXFLAGS ?= -std=c++17 -O3 -march=native FCPW_SIMD_WIDTH ?= 8 INCLUDES := -I$(EBGEOMETRY_HOME) \ - -I$(EBGEOMETRY_HOME)/Submodules/fcpw/include \ - -I$(EBGEOMETRY_HOME)/Submodules/fcpw/deps/eigen \ - -I$(EBGEOMETRY_HOME)/Submodules/fcpw/deps/enoki/include \ - -I$(EBGEOMETRY_HOME)/Submodules/TriangleMeshDistance/TriangleMeshDistance/include + -I$(EBGEOMETRY_HOME)/Benchmark/fcpw/include \ + -I$(EBGEOMETRY_HOME)/Benchmark/fcpw/deps/eigen \ + -I$(EBGEOMETRY_HOME)/Benchmark/fcpw/deps/enoki/include \ + -I$(EBGEOMETRY_HOME)/Benchmark/TriangleMeshDistance/TriangleMeshDistance/include DEFINES := -DFCPW_USE_ENOKI -DFCPW_SIMD_WIDTH=$(FCPW_SIMD_WIDTH) diff --git a/Benchmark/MeshSDF/main.cpp b/Benchmark/MeshSDF/main.cpp index b3dd0c3d..479ea5ef 100644 --- a/Benchmark/MeshSDF/main.cpp +++ b/Benchmark/MeshSDF/main.cpp @@ -42,7 +42,7 @@ using SDF = EBGeometry::TriMeshSDF; int main(int argc, char** argv) { - const std::string objFile = (argc > 1) ? argv[1] : "../../Submodules/common-3d-test-models/data/armadillo.obj"; + const std::string objFile = (argc > 1) ? argv[1] : "../../common-3d-test-models/data/armadillo.obj"; constexpr std::size_t numQueries = 100000; constexpr std::size_t sampleSize = 500; constexpr std::size_t maxLeafGroups = 4; diff --git a/Benchmark/NearestNeighbor/GNUmakefile b/Benchmark/NearestNeighbor/GNUmakefile index b1236284..0007fbe9 100644 --- a/Benchmark/NearestNeighbor/GNUmakefile +++ b/Benchmark/NearestNeighbor/GNUmakefile @@ -1,6 +1,6 @@ # Benchmark: PointCloudBVH vs picoflann vs nanoflann (all-nearest-neighbor). -# Requires the Submodules/{picoflann,nanoflann} submodules: -# git submodule update --init Submodules/picoflann Submodules/nanoflann +# Requires the Benchmark/{picoflann,nanoflann} submodules: +# git submodule update --init Benchmark/picoflann Benchmark/nanoflann # EBGEOMETRY_HOME defaults to the repo root (two levels up). EBGEOMETRY_HOME ?= ../.. @@ -8,8 +8,8 @@ CXX ?= g++ CXXFLAGS ?= -std=c++17 -O3 -march=native INCLUDES := -I$(EBGEOMETRY_HOME) \ - -I$(EBGEOMETRY_HOME)/Submodules/picoflann \ - -I$(EBGEOMETRY_HOME)/Submodules/nanoflann/include + -I$(EBGEOMETRY_HOME)/Benchmark/picoflann \ + -I$(EBGEOMETRY_HOME)/Benchmark/nanoflann/include TARGET := NearestNeighbor.ex diff --git a/Benchmark/README.md b/Benchmark/README.md index 47aa5e0d..d1bdda9b 100644 --- a/Benchmark/README.md +++ b/Benchmark/README.md @@ -1,30 +1,36 @@ Benchmark ========= -Head-to-head performance comparisons of EBGeometry against other open-source geometry-query -libraries. Like `Integrations/`, this is **illustrative and not part of CI** — it depends on -external libraries with their own build systems and release cadences (see issue #109). +Reproducible micro-benchmarks that put EBGeometry's geometry queries next to other open-source +libraries on tasks they have in common, so a reader can see how the tradeoffs play out and re-run +the measurements on their own hardware. The point is the methodology and the cross-checks, not a +scoreboard — numbers vary widely across machines, compilers, and ISAs, so the tables below are +illustrative single-machine snapshots rather than definitive results. -The comparison libraries are pinned as git submodules under [`Submodules/`](../Submodules): +The comparison libraries have their own build systems and release cadences, so they are not linked +into EBGeometry's normal test build. CI *compiles* each benchmark (a smoke test against bit-rot) but +does not *run* or time them (see issue #109). They are pinned as git submodules directly under +`Benchmark/`: * [nanoflann](https://github.com/jlblancoc/nanoflann) — header-only KD-tree (point kNN) * [picoflann](https://github.com/rmsalinas/picoflann) — tiny header-only KD-tree (point kNN) * [fcpw](https://github.com/rohan-sawhney/fcpw) — closest-point / SDF on triangle meshes * [TriangleMeshDistance](https://github.com/InteractiveComputerGraphics/TriangleMeshDistance) — header-only signed distance to triangle meshes -Fetch them (and the mesh submodule) with: +Fetch them (and the top-level mesh submodule) with: ```bash -git submodule update --init Submodules/nanoflann Submodules/picoflann Submodules/fcpw \ - Submodules/TriangleMeshDistance Submodules/common-3d-test-models -git -C Submodules/fcpw submodule update --init deps/eigen # fcpw's Eigen (skip the GPU slang-rhi dep) +git submodule update --init Benchmark/nanoflann Benchmark/picoflann Benchmark/fcpw \ + Benchmark/TriangleMeshDistance common-3d-test-models +git -C Benchmark/fcpw submodule update --init deps/eigen # fcpw's Eigen (skip the GPU slang-rhi dep) ``` fcpw's Enoki CPU-vectorization headers are vendored inside the fcpw submodule, so no extra fetch is needed for it. Each benchmark has a `GNUmakefile` (`make && ./.ex`). Every result is cross-checked against a -brute-force / independent baseline so a wrong answer shows up as a mismatch. +brute-force / independent baseline so a wrong answer shows up as a mismatch — correctness is the part +that transfers across machines even when the timings do not. `NearestNeighbor/` — all-nearest-neighbor on a point cloud ---------------------------------------------------------- @@ -42,11 +48,13 @@ picoflann ~65 0.47 ``` - The KD-tree queries are iterated in **Hilbert order** so their node cache is as warm as - `PointCloudBVH`'s leaf-order batch (natural order is ~2x slower — unfair). That order costs the - KD-trees a one-time ~0.18 us/pt spatial sort; `PointCloudBVH` reuses the ordering its build already - produced, for free. Folding the sort in: `PointCloudBVH` ~0.29 vs nanoflann ~0.43 vs picoflann ~0.65. -- nanoflann has the fastest raw traversal; `PointCloudBVH` wins the end-to-end all-NN job (free query - order + faster build) and builds ~1.5x faster than nanoflann. + `PointCloudBVH`'s leaf-order batch (querying in natural order is ~2x slower, which would not be a + like-for-like comparison). That order costs the KD-trees a one-time ~0.18 us/pt spatial sort; + `PointCloudBVH` reuses the ordering its build already produced. Folding the sort in: + `PointCloudBVH` ~0.29 vs nanoflann ~0.43 vs picoflann ~0.65. +- On this machine nanoflann has the fastest raw per-query traversal, while `PointCloudBVH` comes out + ahead on the end-to-end all-nearest-neighbor job because it gets the query order for free and builds + ~1.5x faster. Which end of that tradeoff matters depends on the workload. `MeshSDF/` — closest-point on a triangle mesh --------------------------------------------- diff --git a/Submodules/TriangleMeshDistance b/Benchmark/TriangleMeshDistance similarity index 100% rename from Submodules/TriangleMeshDistance rename to Benchmark/TriangleMeshDistance diff --git a/Submodules/fcpw b/Benchmark/fcpw similarity index 100% rename from Submodules/fcpw rename to Benchmark/fcpw diff --git a/Submodules/nanoflann b/Benchmark/nanoflann similarity index 100% rename from Submodules/nanoflann rename to Benchmark/nanoflann diff --git a/Submodules/picoflann b/Benchmark/picoflann similarity index 100% rename from Submodules/picoflann rename to Benchmark/picoflann diff --git a/Docs/Sphinx/source/Benchmark.rst b/Docs/Sphinx/source/Benchmark.rst index fed6ae7d..699b1b1e 100644 --- a/Docs/Sphinx/source/Benchmark.rst +++ b/Docs/Sphinx/source/Benchmark.rst @@ -5,23 +5,27 @@ Benchmarks .. important:: - Like the :ref:`Chap:Integrations` examples, the benchmarks are illustrative and **not built or run - as part of EBGeometry's continuous integration** -- they depend on external libraries with their - own build systems. They pin those libraries as git submodules under :file:`Submodules/`. See the - benchmark tracking issue on GitHub for context and planned additions. - -The :file:`Benchmark/` folder compares EBGeometry against other open-source geometry-query libraries -on tasks they have in common. Fetch the comparison libraries (and the mesh submodule) with: + The benchmarks are illustrative, not a scoreboard: reported numbers vary widely across machines, + compilers, and ISAs, so the tables are single-machine snapshots meant to be re-run rather than + quoted. CI *compiles* each benchmark as a smoke test against bit-rot, but does **not** run or time + them -- they depend on external libraries with their own build systems and release cadences. Those + libraries are pinned as git submodules directly under :file:`Benchmark/`. See the benchmark + tracking issue on GitHub for context and planned additions. + +The :file:`Benchmark/` folder places EBGeometry next to other open-source geometry-query libraries on +tasks they have in common, so the tradeoffs are visible and reproducible. Every result is +cross-checked against an independent baseline -- correctness is the part that transfers across +machines even when the timings do not. Fetch the comparison libraries (and the top-level mesh +submodule) with: .. code-block:: bash - git submodule update --init Submodules/nanoflann Submodules/picoflann Submodules/fcpw \ - Submodules/TriangleMeshDistance Submodules/common-3d-test-models - git -C Submodules/fcpw submodule update --init deps/eigen # fcpw's Eigen (skip the GPU dep) + git submodule update --init Benchmark/nanoflann Benchmark/picoflann Benchmark/fcpw \ + Benchmark/TriangleMeshDistance common-3d-test-models + git -C Benchmark/fcpw submodule update --init deps/eigen # fcpw's Eigen (skip the GPU dep) -Each benchmark ships a ``GNUmakefile`` (``make && ./.ex``) and cross-checks every result against -a brute-force / independent baseline. See each folder's ``README.md`` for the full detail and -representative numbers. +Each benchmark ships a ``GNUmakefile`` (``make && ./.ex``). See each folder's ``README.md`` for +the full detail and representative numbers. * :file:`Benchmark/NearestNeighbor` -- all-nearest-neighbor over a point cloud: ``PointCloudBVH`` vs `nanoflann `_ vs diff --git a/Docs/Sphinx/source/ContinuousIntegration.rst b/Docs/Sphinx/source/ContinuousIntegration.rst index 36b93fd7..0cece075 100644 --- a/Docs/Sphinx/source/ContinuousIntegration.rst +++ b/Docs/Sphinx/source/ContinuousIntegration.rst @@ -9,10 +9,11 @@ check: code formatting (``clang-format``) and static analysis (``clang-tidy``, a correctness and assurance (the Catch2 unit-test suite, under multiple compilers, SIMD levels, and both ``float`` and ``double`` precision; every bundled example, built and run via CMake, GNU Make, and direct compiler invocation, under GCC, Clang, and Intel's ``icpx``; AddressSanitizer -and UndefinedBehaviorSanitizer runs of the same test suite); spelling (``codespell``); license -and copyright compliance (REUSE); and the project's documentation (a warnings-as-errors Doxygen -build, and HTML/PDF Sphinx builds). A single aggregator job (``CI-passed``) then gates on all of -the above so branch-protection rules only need to target one required check. +and UndefinedBehaviorSanitizer runs of the same test suite; a compile-only smoke test of the +:file:`Benchmark/` programs against their external comparison libraries); spelling (``codespell``); +license and copyright compliance (REUSE); and the project's documentation (a warnings-as-errors +Doxygen build, and HTML/PDF Sphinx builds). A single aggregator job (``CI-passed``) then gates on +all of the above so branch-protection rules only need to target one required check. .. contents:: On this page :local: @@ -122,6 +123,15 @@ Configures with the ``debug-san`` preset (examples disabled) across a matrix of ``{g++-12, clang++-14}`` × SIMD levels ``{none, avx}``, with ``-DEBGEOMETRY_TEST_BOTH_PRECISIONS=ON``, and runs ``ctest --preset debug-san`` under AddressSanitizer and UndefinedBehaviorSanitizer. +Benchmark-Compile +~~~~~~~~~~~~~~~~~~ + +Fetches each benchmark's comparison-library submodules and *compiles* the two programs under +:file:`Benchmark/` (matrix over ``NearestNeighbor`` and ``MeshSDF``) via their ``GNUmakefile``. This +is a compile-only smoke test guarding against bit-rot as the library evolves -- the benchmarks are +never run or timed in CI, so no mesh data is fetched and no result is treated as a measurement (see +:ref:`Chap:Benchmark`). + CI-passed ~~~~~~~~~ @@ -145,6 +155,7 @@ Dependency graph +-- Unit-Tests +-- Release-Test +-- Sanitizers + +-- Benchmark-Compile (all of the above except Static-analysis) --> CI-passed ``Formatting``, ``Codespell``, ``Reuse``, and ``Doxygen-check`` themselves have no diff --git a/Docs/Sphinx/source/ExampleCSGUnion.rst b/Docs/Sphinx/source/ExampleCSGUnion.rst index 12ab8cc4..483e589b 100644 --- a/Docs/Sphinx/source/ExampleCSGUnion.rst +++ b/Docs/Sphinx/source/ExampleCSGUnion.rst @@ -16,4 +16,4 @@ for how to compile it with CMake, GNU Make, or a direct compiler invocation. cd Examples/CSGUnion ./CSGUnion.ex # defaults to cow.obj - ./CSGUnion.ex ../../Submodules/common-3d-test-models/data/cow.obj + ./CSGUnion.ex ../../common-3d-test-models/data/cow.obj diff --git a/Docs/Sphinx/source/ExampleMeshSDF.rst b/Docs/Sphinx/source/ExampleMeshSDF.rst index 1f0863df..59fb5af4 100644 --- a/Docs/Sphinx/source/ExampleMeshSDF.rst +++ b/Docs/Sphinx/source/ExampleMeshSDF.rst @@ -23,7 +23,7 @@ for how to compile it with CMake, GNU Make, or a direct compiler invocation. cd Examples/MeshSDF ./MeshSDF.ex # defaults to armadillo.obj - ./MeshSDF.ex ../../Submodules/common-3d-test-models/data/cow.obj # or pick another mesh + ./MeshSDF.ex ../../common-3d-test-models/data/cow.obj # or pick another mesh With no argument the example loads ``armadillo.obj`` from the ``common-3d-test-models`` submodule, so make sure it is checked out first (see :ref:`Sec:Cloning`). diff --git a/Docs/Sphinx/source/ObtainingEBGeometry.rst b/Docs/Sphinx/source/ObtainingEBGeometry.rst index 342960b5..10217425 100644 --- a/Docs/Sphinx/source/ObtainingEBGeometry.rst +++ b/Docs/Sphinx/source/ObtainingEBGeometry.rst @@ -13,7 +13,7 @@ Clone the repository from `GitHub `_: The core library is header-only and completely self-contained once cloned. However, the ready-to-run examples in :file:`Examples/` read surface meshes from the `common-3d-test-models `_ collection, -which is bundled as a git submodule (``Submodules/common-3d-test-models/``) at the repository root. +which is bundled as a git submodule (``common-3d-test-models/``) at the repository root. If you intend to run the bundled examples, clone with the submodule in one step instead: .. code-block:: bash @@ -26,13 +26,13 @@ If you already cloned without ``--recurse-submodules``, fetch the submodule afte git submodule update --init --recursive -The meshes are then available as ``.obj`` files under ``Submodules/common-3d-test-models/data/``. +The meshes are then available as ``.obj`` files under ``common-3d-test-models/data/``. Some mesh-based examples take a mesh path on the command line, resolved relative to the run directory (each example is run from its own source folder), for example: .. code-block:: bash - ./a.out ../../Submodules/common-3d-test-models/data/armadillo.obj + ./a.out ../../common-3d-test-models/data/armadillo.obj Running an example with no argument falls back to a default mesh from the submodule, so the submodule must be checked out for the examples to run. diff --git a/Examples/CMakeLists.txt b/Examples/CMakeLists.txt index 86549f74..6defb562 100644 --- a/Examples/CMakeLists.txt +++ b/Examples/CMakeLists.txt @@ -5,7 +5,7 @@ # and layer on the repo-side concerns: optional sanitizers and test registration. # # The binary is run with its source directory as the working directory so that -# relative mesh paths (e.g. "../../Submodules/common-3d-test-models/data/*.obj") +# relative mesh paths (e.g. "../../common-3d-test-models/data/*.obj") # in each example resolve correctly. function(ebgeometry_add_example name) add_subdirectory(${name}) diff --git a/Examples/CSGUnion/README.md b/Examples/CSGUnion/README.md index 05ea1262..c5bcc62f 100644 --- a/Examples/CSGUnion/README.md +++ b/Examples/CSGUnion/README.md @@ -63,7 +63,7 @@ Running Run from this directory so the default mesh path resolves: ./CSGUnion.ex - ./CSGUnion.ex ../../Submodules/common-3d-test-models/data/cow.obj + ./CSGUnion.ex ../../common-3d-test-models/data/cow.obj With no argument the example loads `cow.obj` from the `common-3d-test-models` submodule, so make sure it is checked out diff --git a/Examples/CSGUnion/main.cpp b/Examples/CSGUnion/main.cpp index 15c299db..e7590371 100644 --- a/Examples/CSGUnion/main.cpp +++ b/Examples/CSGUnion/main.cpp @@ -40,7 +40,7 @@ main(int argc, char* argv[]) // Mesh to merge with the sphere. Pass a path on the command line, or fall back to an OBJ file from // the common-3d-test-models submodule (path is relative to this example's source folder, where the // executable is run). See the "Building and using" docs for how to fetch the submodule. - std::string file = "../../Submodules/common-3d-test-models/data/cow.obj"; + std::string file = "../../common-3d-test-models/data/cow.obj"; if (argc >= 2) { file = std::string(argv[1]); } diff --git a/Examples/MeshSDF/README.md b/Examples/MeshSDF/README.md index c60123ab..c81d9b40 100644 --- a/Examples/MeshSDF/README.md +++ b/Examples/MeshSDF/README.md @@ -68,7 +68,7 @@ Running Run from this directory so the default mesh path resolves: ./MeshSDF.ex - ./MeshSDF.ex ../../Submodules/common-3d-test-models/data/cow.obj + ./MeshSDF.ex ../../common-3d-test-models/data/cow.obj With no argument the example loads `armadillo.obj` from the `common-3d-test-models` submodule, so make sure it is checked out diff --git a/Examples/MeshSDF/main.cpp b/Examples/MeshSDF/main.cpp index 1380c2c4..874497b8 100644 --- a/Examples/MeshSDF/main.cpp +++ b/Examples/MeshSDF/main.cpp @@ -31,9 +31,9 @@ int main(int argc, char* argv[]) { // Path to a surface mesh (STL/PLY/VTK/OBJ). Pass one on the command line, e.g. - // ./a.out ../../Submodules/common-3d-test-models/data/cow.obj + // ./a.out ../../common-3d-test-models/data/cow.obj // Paths are resolved relative to the run directory (this example's source folder when run via ctest). - std::string file = "../../Submodules/common-3d-test-models/data/armadillo.obj"; + std::string file = "../../common-3d-test-models/data/armadillo.obj"; if (argc == 2) { file = std::string(argv[1]); diff --git a/Examples/README.md b/Examples/README.md index e7903cfd..bc454aa6 100644 --- a/Examples/README.md +++ b/Examples/README.md @@ -17,5 +17,5 @@ Examples that couple EBGeometry to a third-party application code (AMReX, Chombo [`Integrations/`](../Integrations/README.md) at the repository root instead. The mesh files used by the examples come from the [common-3d-test-models](https://github.com/alecjacobson/common-3d-test-models) -git submodule (`Submodules/common-3d-test-models/`) at the repository root. Clone the repository with `--recurse-submodules` (or run +git submodule (`common-3d-test-models/`) at the repository root. Clone the repository with `--recurse-submodules` (or run `git submodule update --init --recursive`); see the "Building and using" documentation for details. diff --git a/Integrations/AMReX/MeshSDF/main.cpp b/Integrations/AMReX/MeshSDF/main.cpp index 46a89b50..a46146d5 100644 --- a/Integrations/AMReX/MeshSDF/main.cpp +++ b/Integrations/AMReX/MeshSDF/main.cpp @@ -96,7 +96,7 @@ main(int argc, char* argv[]) // Mesh file (STL/PLY/VTK/OBJ). Override with 'filename=' in the inputs file. The default is an OBJ // from the common-3d-test-models submodule (see the "Building and using" docs for how to fetch it); the // path is relative to this example's source folder, where the executable is run. - filename = "../../Submodules/common-3d-test-models/data/armadillo.obj"; + filename = "../../common-3d-test-models/data/armadillo.obj"; pp.query("filename", filename); RealBox rb({-100, -75, -100}, {100, 125, 100}); // bounds the default armadillo mesh diff --git a/Integrations/AMReX/PaintEB/main.cpp b/Integrations/AMReX/PaintEB/main.cpp index 7cf1cf86..bf7cee59 100644 --- a/Integrations/AMReX/PaintEB/main.cpp +++ b/Integrations/AMReX/PaintEB/main.cpp @@ -112,7 +112,7 @@ main(int argc, char* argv[]) // Mesh file (STL/PLY/VTK/OBJ). Override with 'filename=' in the inputs file. The default is a // small OBJ from the common-3d-test-models submodule (see the "Building and using" docs for how to // fetch it); the path is relative to this example's source folder, where the executable is run. - filename = "../../Submodules/common-3d-test-models/data/suzanne.obj"; + filename = "../../common-3d-test-models/data/suzanne.obj"; pp.query("filename", filename); RealBox rb({-5, -1, 2}, {0, 4, 6}); // bounds the default suzanne mesh diff --git a/REUSE.toml b/REUSE.toml index 08c965ab..1f250d40 100644 --- a/REUSE.toml +++ b/REUSE.toml @@ -6,7 +6,7 @@ # Source code carries inline SPDX headers instead; see Source/ and Examples/*/main.cpp. # # NOTE: third-party example meshes are no longer redistributed here; they are pulled in -# via the common-3d-test-models git submodule under Submodules/ (a separate repository +# via the top-level common-3d-test-models git submodule (a separate repository # with its own licensing, which REUSE does not descend into). version = 1 diff --git a/Submodules/common-3d-test-models b/common-3d-test-models similarity index 100% rename from Submodules/common-3d-test-models rename to common-3d-test-models From 2ddfc6ce76d25000e0f4f51ece035356bc4a8307 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 15 Jul 2026 21:35:16 +0000 Subject: [PATCH 05/17] Add PointCloudHashGrid and a sphere-surface distribution to the NN benchmark Extend Benchmark/NearestNeighbor to compare EBGeometry's PointCloudHashGrid alongside PointCloudBVH, picoflann, and nanoflann, and run the whole comparison over two point distributions in the same program: the existing uniform unit-cube cloud and a new uniformly-sampled unit-sphere surface (a 2D manifold that is locally dense but globally hollow -- a harder case for a uniform grid). Every method is still cross-checked against the brute-force baseline on both cases. Update Benchmark/README.md's NearestNeighbor section to describe the two distributions and the four methods, and to record the build/query tradeoff the hash grid exposes. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01JhdRJydutVJqQjg2Fuee9Z --- Benchmark/NearestNeighbor/main.cpp | 113 ++++++++++++++++++++++------- Benchmark/README.md | 41 +++++++---- 2 files changed, 114 insertions(+), 40 deletions(-) diff --git a/Benchmark/NearestNeighbor/main.cpp b/Benchmark/NearestNeighbor/main.cpp index 9b1bcb50..69bc1153 100644 --- a/Benchmark/NearestNeighbor/main.cpp +++ b/Benchmark/NearestNeighbor/main.cpp @@ -2,9 +2,15 @@ // // SPDX-License-Identifier: GPL-3.0-or-later -// Benchmark: EBGeometry PointCloudBVH vs picoflann vs nanoflann, all-nearest-neighbor. -// 500k random 3D points in the unit cube (double). Every point's nearest OTHER point. All three +// Benchmark: EBGeometry PointCloudBVH and PointCloudHashGrid vs picoflann vs nanoflann, +// all-nearest-neighbor. 500k 3D points (double). Every point's nearest OTHER point. All four // verified against a brute-force sample. Both KD-trees used vanilla; distances are squared throughout. +// +// Two point distributions are benchmarked in turn, since spatial data structures behave very +// differently depending on how the points fill space: +// 1. Uniform in the unit cube -- points fill a 3D volume evenly (the easy, balanced case). +// 2. On the unit-sphere surface -- points lie on a 2D manifold: locally dense, globally hollow, +// a harder case for uniform grids (many empty cells inside). #include #include @@ -13,6 +19,8 @@ #include #include #include +#include +#include #include #include @@ -62,6 +70,31 @@ struct NanoCloud }; using NanoTree = nanoflann::KDTreeSingleIndexAdaptor, NanoCloud, 3>; +// Uniformly sample the surface of the unit sphere: draw a 3D Gaussian and normalize, which is +// rotationally symmetric and hence uniform over the sphere. Points lie on a 2D manifold embedded in +// 3D -- locally dense, globally hollow. +std::vector +samplePointsOnSphere(std::size_t a_count, std::uint64_t a_seed) +{ + std::mt19937_64 rng(a_seed); + std::normal_distribution gauss(T(0), T(1)); + + std::vector points; + points.reserve(a_count); + + for (std::size_t i = 0; i < a_count; i++) { + Vec3 v(gauss(rng), gauss(rng), gauss(rng)); + T len = v.length(); + if (len < std::numeric_limits::min()) { + v = Vec3(T(1), T(0), T(0)); // Degenerate zero draw; nudge onto the sphere. + len = T(1); + } + points.emplace_back(v / len); + } + + return points; +} + T bruteForceNN2(std::size_t a_self, const std::vector& a_pos) { @@ -74,32 +107,30 @@ bruteForceNN2(std::size_t a_self, const std::vector& a_pos) return best; } -} // namespace - -int -main() +// Run the full all-nearest-neighbor comparison over one point distribution and print a table. +void +runCase(const std::string& a_label, const std::vector& a_positions) { - std::cout << "All-nearest-neighbor: PointCloudBVH vs picoflann vs nanoflann\n"; - std::cout << " Points = " << numPoints << " (double, unit cube)\n\n"; - - const std::vector positions = EBGeometry::Random::samplePoints(numPoints, pointSeed); - const std::vector meta(numPoints); + const std::size_t n = a_positions.size(); + const std::vector meta(n); EBGeometry::SimpleTimer timer; + std::cout << "== " << a_label << " (" << n << " points, double) ==\n"; + // Query the flann trees in a spatially-coherent (Hilbert) order too, so their node cache is as warm // as EBGeometry's leaf-order batch. EBGeometry gets its order free from the build; the flann libs // must sort -- time that once so it can be folded in if desired. timer.start(); - const std::vector order = EBGeometry::SFC::order(positions); + const std::vector order = EBGeometry::SFC::order(a_positions); timer.stop(); - const double sortUsPerPt = 1.0e6 * timer.seconds() / double(numPoints); + const double sortUsPerPt = 1.0e6 * timer.seconds() / double(n); - const std::size_t stride = numPoints / sampleSize; + const std::size_t stride = n / sampleSize; std::vector truth(sampleSize); timer.start(); for (std::size_t s = 0; s < sampleSize; s++) { - truth[s] = bruteForceNN2(s * stride, positions); + truth[s] = bruteForceNN2(s * stride, a_positions); } timer.stop(); const double bruteUsPerPt = 1.0e6 * timer.seconds() / double(sampleSize); @@ -124,14 +155,14 @@ main() // ── EBGeometry PointCloudBVH (batched all-NN) ── { timer.start(); - const EBGeometry::PointCloudBVH bvh(positions, meta); + const EBGeometry::PointCloudBVH bvh(a_positions, meta); timer.stop(); const double buildMs = 1.0e3 * timer.seconds(); timer.start(); const auto graph = bvh.allNearestNeighbors(1); timer.stop(); - const double queryUs = 1.0e6 * timer.seconds() / double(numPoints); + const double queryUs = 1.0e6 * timer.seconds() / double(n); std::size_t bad = 0; for (std::size_t s = 0; s < sampleSize; s++) { @@ -140,18 +171,37 @@ main() row("PointCloudBVH", buildMs, queryUs, bad); } + // ── EBGeometry PointCloudHashGrid (batched all-NN, uniform grid) ── + { + timer.start(); + const EBGeometry::PointCloudHashGrid grid(a_positions, meta); + timer.stop(); + const double buildMs = 1.0e3 * timer.seconds(); + + timer.start(); + const auto graph = grid.allNearestNeighbors(1); + timer.stop(); + const double queryUs = 1.0e6 * timer.seconds() / double(n); + + std::size_t bad = 0; + for (std::size_t s = 0; s < sampleSize; s++) { + bad += !ok(graph[s * stride].distanceSquared, s); + } + row("PointCloudHashGrid", buildMs, queryUs, bad); + } + // ── picoflann (per-point searchKnn) ── { picoflann::KdTreeIndex<3, Vec3Adapter> kdtree; timer.start(); - kdtree.build(positions); + kdtree.build(a_positions); timer.stop(); const double buildMs = 1.0e3 * timer.seconds(); volatile T sink = T(0); timer.start(); for (const std::uint32_t p : order) { - const auto res = kdtree.searchKnn(positions, positions[p], 2); + const auto res = kdtree.searchKnn(a_positions, a_positions[p], 2); for (const auto& pr : res) { if (pr.first != p) { sink += pr.second; @@ -161,11 +211,11 @@ main() } timer.stop(); (void)sink; - const double queryUs = 1.0e6 * timer.seconds() / double(numPoints); + const double queryUs = 1.0e6 * timer.seconds() / double(n); std::size_t bad = 0; for (std::size_t s = 0; s < sampleSize; s++) { - const auto res = kdtree.searchKnn(positions, positions[s * stride], 2); + const auto res = kdtree.searchKnn(a_positions, a_positions[s * stride], 2); T got = std::numeric_limits::max(); for (const auto& pr : res) { if (pr.first != s * stride) { @@ -180,7 +230,7 @@ main() // ── nanoflann (per-point findNeighbors, k=2) ── { - NanoCloud cloud{positions}; + NanoCloud cloud{a_positions}; timer.start(); NanoTree index(3, cloud, nanoflann::KDTreeSingleIndexAdaptorParams(10 /* leaf_max_size */)); index.buildIndex(); @@ -192,7 +242,7 @@ main() T d2[2]; nanoflann::KNNResultSet rs(2); rs.init(idx, d2); - const T qp[3] = {positions[i][0], positions[i][1], positions[i][2]}; + const T qp[3] = {a_positions[i][0], a_positions[i][1], a_positions[i][2]}; index.findNeighbors(rs, qp); // idx[0] is the point itself (distance 0); take the first neighbor that is not itself. return (idx[0] != i) ? d2[0] : d2[1]; @@ -205,7 +255,7 @@ main() } timer.stop(); (void)sink; - const double queryUs = 1.0e6 * timer.seconds() / double(numPoints); + const double queryUs = 1.0e6 * timer.seconds() / double(n); std::size_t bad = 0; for (std::size_t s = 0; s < sampleSize; s++) { @@ -214,9 +264,20 @@ main() row("nanoflann", buildMs, queryUs, bad); } - std::cout << "\n flann query loops iterate in Hilbert order (warm cache, like EBGeometry's leaf order).\n"; + std::cout << " flann query loops iterate in Hilbert order (warm cache, like EBGeometry's leaf order).\n"; std::cout << " One-time Hilbert sort the flann libs need for that order: " << std::setprecision(3) << sortUsPerPt - << " us/pt (add to their query if counted;\n EBGeometry reuses its build order for free).\n"; + << " us/pt\n (add to their query if counted; EBGeometry reuses its build order for free).\n\n"; +} + +} // namespace + +int +main() +{ + std::cout << "All-nearest-neighbor: PointCloudBVH & PointCloudHashGrid vs picoflann vs nanoflann\n\n"; + + runCase("Uniform in the unit cube", EBGeometry::Random::samplePoints(numPoints, pointSeed)); + runCase("On the unit-sphere surface", samplePointsOnSphere(numPoints, pointSeed)); return 0; } diff --git a/Benchmark/README.md b/Benchmark/README.md index d1bdda9b..43cb752d 100644 --- a/Benchmark/README.md +++ b/Benchmark/README.md @@ -35,26 +35,39 @@ that transfers across machines even when the timings do not. `NearestNeighbor/` — all-nearest-neighbor on a point cloud ---------------------------------------------------------- -For every point in a 500,000-point uniform cloud in the unit cube (double precision), find its -nearest *other* point. Compares `PointCloudBVH` vs picoflann vs nanoflann. +For every point in a 500,000-point cloud (double precision), find its nearest *other* point. Compares +EBGeometry's `PointCloudBVH` and `PointCloudHashGrid` against picoflann and nanoflann. The same +comparison is run over two point distributions, since spatial structures behave very differently +depending on how the points fill space: -Representative result (one machine): +* **Uniform in the unit cube** — points fill a 3D volume evenly (the balanced, easy case). +* **On the unit-sphere surface** — points lie on a 2D manifold: locally dense, globally hollow. This + is the harder case for a uniform grid, whose bounding box is then mostly empty interior cells. + +Representative result (one machine, illustrative — see the note on machine dependence above): ``` -Method Build(ms) Query(us/pt) -PointCloudBVH ~78 0.29 -nanoflann ~120 0.25 -picoflann ~65 0.47 + uniform cube sphere surface +Method Build(ms) Query(us/pt) Build(ms) Query(us/pt) +PointCloudBVH ~165 0.90 ~155 0.66 +PointCloudHashGrid ~19 1.85 ~17 2.04 +picoflann ~132 0.87 ~130 0.55 +nanoflann ~262 0.44 ~273 0.33 ``` - The KD-tree queries are iterated in **Hilbert order** so their node cache is as warm as - `PointCloudBVH`'s leaf-order batch (querying in natural order is ~2x slower, which would not be a - like-for-like comparison). That order costs the KD-trees a one-time ~0.18 us/pt spatial sort; - `PointCloudBVH` reuses the ordering its build already produced. Folding the sort in: - `PointCloudBVH` ~0.29 vs nanoflann ~0.43 vs picoflann ~0.65. -- On this machine nanoflann has the fastest raw per-query traversal, while `PointCloudBVH` comes out - ahead on the end-to-end all-nearest-neighbor job because it gets the query order for free and builds - ~1.5x faster. Which end of that tradeoff matters depends on the workload. + EBGeometry's leaf-order batch (querying in natural order is ~2x slower, which would not be a + like-for-like comparison). That order costs the KD-trees a one-time spatial sort (~0.35 us/pt here); + the EBGeometry structures reuse the ordering their build already produced. +- **`PointCloudHashGrid` trades query speed for build speed**: an O(N) uniform grid builds ~8x faster + than the BVH but scans neighbor cells per query, so it queries ~2x slower. It is the weakest on the + sphere surface (query ~2.0 us/pt) — the hollow distribution leaves its grid mostly empty while the + occupied surface cells are denser than the ~1-point-per-cell target. +- The tree/BVH methods, by contrast, get *faster* on the sphere surface than in the cube (the local + neighborhood is effectively lower-dimensional, so pruning is tighter). nanoflann has the fastest + raw per-query traversal throughout; `PointCloudBVH` is competitive end-to-end because it gets its + query order for free and builds faster. Which structure to pick depends on the build/query balance + and the point distribution — that is the point of running both cases. `MeshSDF/` — closest-point on a triangle mesh --------------------------------------------- From dc03101321621029d0c034924d61d88c6e5ae142 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 15 Jul 2026 22:10:57 +0000 Subject: [PATCH 06/17] Prototype eps-approximate pruning in PointCloudBVH and benchmark it Following the nanoflann traversal study: nanoflann prunes the far child with a (1+eps)-loosened bound for approximate search. Add the equivalent to PointCloudBVH as an optional, backward-compatible knob -- allNearestNeighbors gains an a_eps parameter (default 0 = exact, unchanged behavior) that is converted to a squared-distance prune-scale 1/(1+eps)^2 and threaded through query() into both the seeded scalar DFS and the pruneTraverse-based paths. A returned neighbor is then guaranteed within a factor (1+eps) of the true nearest; best.distanceSquared always remains the true distance to a real point. Add a correctness test asserting the (1+eps) guarantee (and that eps==0 reproduces the exact search), and an eps-sweep section to the NearestNeighbor benchmark reporting query time, speedup, %exact, and worst-case distance ratio on both point distributions. Finding: the speedup is modest (~1.1-1.35x even at eps=2) because the all-NN path already seeds each query's bound from the point's own leaf, so the bound starts near-tight and there is little marginal pruning left for eps to remove -- unlike the loose-initial-bound high-dimensional kNN case where eps pays off. The separately-considered no-sort idea needs no prototype: the seeded all-NN path already uses an unordered scalar DFS, so the sort in pruneTraverse never runs on this workload. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01JhdRJydutVJqQjg2Fuee9Z --- Benchmark/NearestNeighbor/main.cpp | 48 ++++++++++++++++++++++- Source/EBGeometry_PointCloudBVH.hpp | 13 ++++-- Source/EBGeometry_PointCloudBVHImplem.hpp | 33 ++++++++++++---- Tests/TestPointCloudBVH.cpp | 25 ++++++++++++ 4 files changed, 107 insertions(+), 12 deletions(-) diff --git a/Benchmark/NearestNeighbor/main.cpp b/Benchmark/NearestNeighbor/main.cpp index 69bc1153..7b14a2df 100644 --- a/Benchmark/NearestNeighbor/main.cpp +++ b/Benchmark/NearestNeighbor/main.cpp @@ -266,7 +266,53 @@ runCase(const std::string& a_label, const std::vector& a_positions) std::cout << " flann query loops iterate in Hilbert order (warm cache, like EBGeometry's leaf order).\n"; std::cout << " One-time Hilbert sort the flann libs need for that order: " << std::setprecision(3) << sortUsPerPt - << " us/pt\n (add to their query if counted; EBGeometry reuses its build order for free).\n\n"; + << " us/pt\n (add to their query if counted; EBGeometry reuses its build order for free).\n"; + + // ── PointCloudBVH eps-approximate all-NN: speed vs measured accuracy ── + // The nanoflann-inspired experiment: loosen the traversal's prune bound by a factor (1+eps) so + // boxes that could only marginally improve the result are skipped. eps == 0 is the exact search + // (identical to the table above). Accuracy is measured against the brute-force truth[] sample: + // "%exact" is how often the returned neighbor is still the true nearest; "max ratio" is the worst + // returned/true distance ratio over the sample (guaranteed <= 1+eps). + { + const EBGeometry::PointCloudBVH bvh(a_positions, meta); + + std::cout << "\n PointCloudBVH eps-approximate all-NN (nanoflann-style loosened prune bound):\n"; + std::cout << " " << std::left << std::setw(8) << "eps" << std::right << std::setw(14) << "Query(us/pt)" + << std::setw(10) << "speedup" << std::setw(9) << "%exact" << std::setw(12) << "max ratio" << '\n'; + + double exactQueryUs = 1.0; + for (const double eps : {0.0, 0.1, 0.25, 0.5, 1.0, 2.0}) { + timer.start(); + const auto graph = bvh.allNearestNeighbors(1, T(eps)); + timer.stop(); + const double queryUs = 1.0e6 * timer.seconds() / double(n); + + if (eps == 0.0) { + exactQueryUs = queryUs; + } + + std::size_t exact = 0; + double maxRatio = 1.0; + for (std::size_t s = 0; s < sampleSize; s++) { + const T got = graph[s * stride].distanceSquared; + + if (ok(got, s)) { + exact++; + } + else { + maxRatio = std::max(maxRatio, std::sqrt(double(got) / std::max(double(truth[s]), 1.0e-300))); + } + } + + std::cout << " " << std::left << std::setw(8) << std::setprecision(2) << eps << std::right << std::setw(14) + << std::setprecision(3) << queryUs << std::setw(9) << std::setprecision(2) << exactQueryUs / queryUs + << "x" << std::setw(8) << std::setprecision(1) << 100.0 * double(exact) / double(sampleSize) << "%" + << std::setw(12) << std::setprecision(4) << maxRatio << '\n'; + } + } + + std::cout << '\n'; } } // namespace diff --git a/Source/EBGeometry_PointCloudBVH.hpp b/Source/EBGeometry_PointCloudBVH.hpp index fd7c9b19..707da9bd 100644 --- a/Source/EBGeometry_PointCloudBVH.hpp +++ b/Source/EBGeometry_PointCloudBVH.hpp @@ -154,11 +154,14 @@ class PointCloudBVH * its own leaf, so the whole batch is cheaper than the sum of independent queries. Result is * flattened row-major: entry [i*a_k + j] is the j-th nearest neighbor of point i (ascending by * distance). - * @param[in] a_k Number of neighbors per point. + * @param[in] a_k Number of neighbors per point. + * @param[in] a_eps Approximation tolerance (>= 0). 0 (the default) is an exact search; a positive + * value returns, for each point, neighbors within a factor (1+a_eps) of the true + * nearest distance, pruning more aggressively for a faster query. * @return A vector of size numPoints()*a_k of Hits. */ [[nodiscard]] inline std::vector - allNearestNeighbors(std::size_t a_k = 1) const; + allNearestNeighbors(std::size_t a_k = 1, T a_eps = T(0)) const; /** * @brief Brute-force closest point to an arbitrary query point (O(N) reference for closestPoint()). @@ -285,6 +288,9 @@ class PointCloudBVH * @param[in] a_seedOff Group offset of the own leaf to seed from (self-queries; ignored if * a_seedCnt == 0). * @param[in] a_seedCnt Group count of the own leaf; 0 means "no seed" (external queries). + * @param[in] a_pruneScale Multiplicative acceptance-bound factor in (0, 1]; 1 (the default) is an + * exact query. A value 1/(1+eps)^2 yields an eps-approximate query -- boxes + * that could improve the result by less than a factor (1+eps) are pruned. */ inline void query(const Vec3T& a_query, @@ -293,7 +299,8 @@ class PointCloudBVH std::size_t& a_found, std::size_t a_exclude, std::uint32_t a_seedOff, - std::uint32_t a_seedCnt) const noexcept; + std::uint32_t a_seedCnt, + T a_pruneScale = T(1)) const noexcept; /** * @brief Brute-force single nearest by full scan (shared by closestPointBruteForce / diff --git a/Source/EBGeometry_PointCloudBVHImplem.hpp b/Source/EBGeometry_PointCloudBVHImplem.hpp index 25d6213b..22cf9c27 100644 --- a/Source/EBGeometry_PointCloudBVHImplem.hpp +++ b/Source/EBGeometry_PointCloudBVHImplem.hpp @@ -240,7 +240,8 @@ PointCloudBVH::query(const Vec3T& a_query, std::size_t& a_found, std::size_t a_exclude, std::uint32_t a_seedOff, - std::uint32_t a_seedCnt) const noexcept + std::uint32_t a_seedCnt, + T a_pruneScale) const noexcept { EBGEOMETRY_EXPECT(a_k >= 1); EBGEOMETRY_EXPECT(a_out != nullptr); @@ -308,8 +309,14 @@ PointCloudBVH::query(const Vec3T& a_query, while (stackTop > 0) { const Node& node = this->m_linearNodes[stack[--stackTop]]; - if (node.getDistanceToBoundingVolume2(a_query) >= best.distanceSquared) { - continue; // stale: best tightened since this node was pushed + // a_pruneScale (== 1/(1+eps)^2, defaulting to 1 for an exact query) shrinks the acceptance + // bound: a box is explored only if it could hold a point closer than best/(1+eps), so an + // eps-approximate query skips boxes whose contents could improve the result by less than that + // factor. best.distanceSquared always stays the true squared distance to a real point. + const T pruneBound = best.distanceSquared * a_pruneScale; + + if (node.getDistanceToBoundingVolume2(a_query) >= pruneBound) { + continue; // stale: best tightened since this node was pushed, or eps-pruned } if (node.isLeaf()) { @@ -323,7 +330,7 @@ PointCloudBVH::query(const Vec3T& a_query, const auto& childOffsets = node.getChildOffsets(); for (std::size_t k = 0; k < K; k++) { - if (this->m_linearNodes[childOffsets[k]].getDistanceToBoundingVolume2(a_query) < best.distanceSquared) { + if (this->m_linearNodes[childOffsets[k]].getDistanceToBoundingVolume2(a_query) < pruneBound) { EBGEOMETRY_EXPECT(stackTop < maxStack); stack[stackTop++] = childOffsets[k]; @@ -342,7 +349,9 @@ PointCloudBVH::query(const Vec3T& a_query, scanLeafBest(a_best, a_off, a_cnt); }; - const auto pruneDist2 = [](const Best& a_best) noexcept -> T { return a_best.distanceSquared; }; + const auto pruneDist2 = [a_pruneScale](const Best& a_best) noexcept -> T { + return a_best.distanceSquared * a_pruneScale; + }; this->pruneTraverse(a_query, best, evalLeaf, pruneDist2); } @@ -419,7 +428,9 @@ PointCloudBVH::query(const Vec3T& a_query, processLeaf(a_state, a_off, a_cnt); }; - const auto pruneDist2 = [](const QState& a_state) noexcept -> T { return a_state.bound; }; + const auto pruneDist2 = [a_pruneScale](const QState& a_state) noexcept -> T { + return a_state.bound * a_pruneScale; + }; this->pruneTraverse(a_query, state, evalLeaf, pruneDist2); @@ -483,12 +494,18 @@ PointCloudBVH::nearestNeighbors(std::size_t a_point, std::size_t template inline std::vector::Hit> -PointCloudBVH::allNearestNeighbors(std::size_t a_k) const +PointCloudBVH::allNearestNeighbors(std::size_t a_k, T a_eps) const { EBGEOMETRY_EXPECT(a_k >= 1); + EBGEOMETRY_EXPECT(a_eps >= T(0)); const std::size_t numPoints = m_positions.size(); + // Approximate-query pruning factor: a box is explored only if it could beat best/(1+eps), so any + // returned neighbor is within a factor (1+a_eps) of the true nearest. a_eps == 0 (the default) makes + // the factor 1 and the search exact. Distances are squared throughout, hence the squared reciprocal. + const T pruneScale = T(1) / ((T(1) + a_eps) * (T(1) + a_eps)); + std::vector result(numPoints * a_k); // Process points in leaf (build) order: the top-down build already grouped them spatially, so @@ -500,7 +517,7 @@ PointCloudBVH::allNearestNeighbors(std::size_t a_k) const std::size_t found = 0; - this->query(m_positions[p], a_k, &result[p * a_k], found, p, m_leafOff[p], m_leafCnt[p]); + this->query(m_positions[p], a_k, &result[p * a_k], found, p, m_leafOff[p], m_leafCnt[p], pruneScale); } return result; diff --git a/Tests/TestPointCloudBVH.cpp b/Tests/TestPointCloudBVH.cpp index 21d4a6e9..7b1c43ac 100644 --- a/Tests/TestPointCloudBVH.cpp +++ b/Tests/TestPointCloudBVH.cpp @@ -146,6 +146,31 @@ TEMPLATE_TEST_CASE("PointCloudBVH queries match brute force", "[PointCloudBVH]", } } + SECTION("eps-approximate allNearestNeighbors stays within the (1+eps) distance guarantee") + { + const auto exact = bvh.allNearestNeighbors(1); + REQUIRE(exact.size() == n); + + for (const T eps : {T(0.25), T(1.0), T(4.0)}) { + const T ratioBound = (T(1) + eps) * (T(1) + eps); // guarantee is on distance, we compare squares + const auto approx = bvh.allNearestNeighbors(1, eps); + REQUIRE(approx.size() == n); + + for (std::size_t i = 0; i < n; i++) { + // The returned neighbor is a real point, so it is never closer than the true nearest, and no + // farther than a factor (1+eps): exactTrue <= approx <= (1+eps)^2 * exactTrue (squared). + CHECK(approx[i].distanceSquared >= exact[i].distanceSquared * (T(1) - looseMargin())); + CHECK(approx[i].distanceSquared <= exact[i].distanceSquared * ratioBound * (T(1) + looseMargin())); + } + } + + // eps == 0 must reproduce the exact search bit-for-bit (the default path). + const auto exact0 = bvh.allNearestNeighbors(1, T(0)); + for (std::size_t i = 0; i < n; i++) { + CHECK_THAT(exact0[i].distanceSquared, withinAbsT(exact[i].distanceSquared, tol)); + } + } + SECTION("accessors return the stored cloud data") { for (std::size_t i = 0; i < n; i += 313) { From 5ec6e96b33ff81e55e5518211b9bec82bcb78af5 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 15 Jul 2026 22:13:14 +0000 Subject: [PATCH 07/17] Revert "Prototype eps-approximate pruning in PointCloudBVH and benchmark it" This reverts commit dc03101321621029d0c034924d61d88c6e5ae142. --- Benchmark/NearestNeighbor/main.cpp | 48 +---------------------- Source/EBGeometry_PointCloudBVH.hpp | 13 ++---- Source/EBGeometry_PointCloudBVHImplem.hpp | 33 ++++------------ Tests/TestPointCloudBVH.cpp | 25 ------------ 4 files changed, 12 insertions(+), 107 deletions(-) diff --git a/Benchmark/NearestNeighbor/main.cpp b/Benchmark/NearestNeighbor/main.cpp index 7b14a2df..69bc1153 100644 --- a/Benchmark/NearestNeighbor/main.cpp +++ b/Benchmark/NearestNeighbor/main.cpp @@ -266,53 +266,7 @@ runCase(const std::string& a_label, const std::vector& a_positions) std::cout << " flann query loops iterate in Hilbert order (warm cache, like EBGeometry's leaf order).\n"; std::cout << " One-time Hilbert sort the flann libs need for that order: " << std::setprecision(3) << sortUsPerPt - << " us/pt\n (add to their query if counted; EBGeometry reuses its build order for free).\n"; - - // ── PointCloudBVH eps-approximate all-NN: speed vs measured accuracy ── - // The nanoflann-inspired experiment: loosen the traversal's prune bound by a factor (1+eps) so - // boxes that could only marginally improve the result are skipped. eps == 0 is the exact search - // (identical to the table above). Accuracy is measured against the brute-force truth[] sample: - // "%exact" is how often the returned neighbor is still the true nearest; "max ratio" is the worst - // returned/true distance ratio over the sample (guaranteed <= 1+eps). - { - const EBGeometry::PointCloudBVH bvh(a_positions, meta); - - std::cout << "\n PointCloudBVH eps-approximate all-NN (nanoflann-style loosened prune bound):\n"; - std::cout << " " << std::left << std::setw(8) << "eps" << std::right << std::setw(14) << "Query(us/pt)" - << std::setw(10) << "speedup" << std::setw(9) << "%exact" << std::setw(12) << "max ratio" << '\n'; - - double exactQueryUs = 1.0; - for (const double eps : {0.0, 0.1, 0.25, 0.5, 1.0, 2.0}) { - timer.start(); - const auto graph = bvh.allNearestNeighbors(1, T(eps)); - timer.stop(); - const double queryUs = 1.0e6 * timer.seconds() / double(n); - - if (eps == 0.0) { - exactQueryUs = queryUs; - } - - std::size_t exact = 0; - double maxRatio = 1.0; - for (std::size_t s = 0; s < sampleSize; s++) { - const T got = graph[s * stride].distanceSquared; - - if (ok(got, s)) { - exact++; - } - else { - maxRatio = std::max(maxRatio, std::sqrt(double(got) / std::max(double(truth[s]), 1.0e-300))); - } - } - - std::cout << " " << std::left << std::setw(8) << std::setprecision(2) << eps << std::right << std::setw(14) - << std::setprecision(3) << queryUs << std::setw(9) << std::setprecision(2) << exactQueryUs / queryUs - << "x" << std::setw(8) << std::setprecision(1) << 100.0 * double(exact) / double(sampleSize) << "%" - << std::setw(12) << std::setprecision(4) << maxRatio << '\n'; - } - } - - std::cout << '\n'; + << " us/pt\n (add to their query if counted; EBGeometry reuses its build order for free).\n\n"; } } // namespace diff --git a/Source/EBGeometry_PointCloudBVH.hpp b/Source/EBGeometry_PointCloudBVH.hpp index 707da9bd..fd7c9b19 100644 --- a/Source/EBGeometry_PointCloudBVH.hpp +++ b/Source/EBGeometry_PointCloudBVH.hpp @@ -154,14 +154,11 @@ class PointCloudBVH * its own leaf, so the whole batch is cheaper than the sum of independent queries. Result is * flattened row-major: entry [i*a_k + j] is the j-th nearest neighbor of point i (ascending by * distance). - * @param[in] a_k Number of neighbors per point. - * @param[in] a_eps Approximation tolerance (>= 0). 0 (the default) is an exact search; a positive - * value returns, for each point, neighbors within a factor (1+a_eps) of the true - * nearest distance, pruning more aggressively for a faster query. + * @param[in] a_k Number of neighbors per point. * @return A vector of size numPoints()*a_k of Hits. */ [[nodiscard]] inline std::vector - allNearestNeighbors(std::size_t a_k = 1, T a_eps = T(0)) const; + allNearestNeighbors(std::size_t a_k = 1) const; /** * @brief Brute-force closest point to an arbitrary query point (O(N) reference for closestPoint()). @@ -288,9 +285,6 @@ class PointCloudBVH * @param[in] a_seedOff Group offset of the own leaf to seed from (self-queries; ignored if * a_seedCnt == 0). * @param[in] a_seedCnt Group count of the own leaf; 0 means "no seed" (external queries). - * @param[in] a_pruneScale Multiplicative acceptance-bound factor in (0, 1]; 1 (the default) is an - * exact query. A value 1/(1+eps)^2 yields an eps-approximate query -- boxes - * that could improve the result by less than a factor (1+eps) are pruned. */ inline void query(const Vec3T& a_query, @@ -299,8 +293,7 @@ class PointCloudBVH std::size_t& a_found, std::size_t a_exclude, std::uint32_t a_seedOff, - std::uint32_t a_seedCnt, - T a_pruneScale = T(1)) const noexcept; + std::uint32_t a_seedCnt) const noexcept; /** * @brief Brute-force single nearest by full scan (shared by closestPointBruteForce / diff --git a/Source/EBGeometry_PointCloudBVHImplem.hpp b/Source/EBGeometry_PointCloudBVHImplem.hpp index 22cf9c27..25d6213b 100644 --- a/Source/EBGeometry_PointCloudBVHImplem.hpp +++ b/Source/EBGeometry_PointCloudBVHImplem.hpp @@ -240,8 +240,7 @@ PointCloudBVH::query(const Vec3T& a_query, std::size_t& a_found, std::size_t a_exclude, std::uint32_t a_seedOff, - std::uint32_t a_seedCnt, - T a_pruneScale) const noexcept + std::uint32_t a_seedCnt) const noexcept { EBGEOMETRY_EXPECT(a_k >= 1); EBGEOMETRY_EXPECT(a_out != nullptr); @@ -309,14 +308,8 @@ PointCloudBVH::query(const Vec3T& a_query, while (stackTop > 0) { const Node& node = this->m_linearNodes[stack[--stackTop]]; - // a_pruneScale (== 1/(1+eps)^2, defaulting to 1 for an exact query) shrinks the acceptance - // bound: a box is explored only if it could hold a point closer than best/(1+eps), so an - // eps-approximate query skips boxes whose contents could improve the result by less than that - // factor. best.distanceSquared always stays the true squared distance to a real point. - const T pruneBound = best.distanceSquared * a_pruneScale; - - if (node.getDistanceToBoundingVolume2(a_query) >= pruneBound) { - continue; // stale: best tightened since this node was pushed, or eps-pruned + if (node.getDistanceToBoundingVolume2(a_query) >= best.distanceSquared) { + continue; // stale: best tightened since this node was pushed } if (node.isLeaf()) { @@ -330,7 +323,7 @@ PointCloudBVH::query(const Vec3T& a_query, const auto& childOffsets = node.getChildOffsets(); for (std::size_t k = 0; k < K; k++) { - if (this->m_linearNodes[childOffsets[k]].getDistanceToBoundingVolume2(a_query) < pruneBound) { + if (this->m_linearNodes[childOffsets[k]].getDistanceToBoundingVolume2(a_query) < best.distanceSquared) { EBGEOMETRY_EXPECT(stackTop < maxStack); stack[stackTop++] = childOffsets[k]; @@ -349,9 +342,7 @@ PointCloudBVH::query(const Vec3T& a_query, scanLeafBest(a_best, a_off, a_cnt); }; - const auto pruneDist2 = [a_pruneScale](const Best& a_best) noexcept -> T { - return a_best.distanceSquared * a_pruneScale; - }; + const auto pruneDist2 = [](const Best& a_best) noexcept -> T { return a_best.distanceSquared; }; this->pruneTraverse(a_query, best, evalLeaf, pruneDist2); } @@ -428,9 +419,7 @@ PointCloudBVH::query(const Vec3T& a_query, processLeaf(a_state, a_off, a_cnt); }; - const auto pruneDist2 = [a_pruneScale](const QState& a_state) noexcept -> T { - return a_state.bound * a_pruneScale; - }; + const auto pruneDist2 = [](const QState& a_state) noexcept -> T { return a_state.bound; }; this->pruneTraverse(a_query, state, evalLeaf, pruneDist2); @@ -494,18 +483,12 @@ PointCloudBVH::nearestNeighbors(std::size_t a_point, std::size_t template inline std::vector::Hit> -PointCloudBVH::allNearestNeighbors(std::size_t a_k, T a_eps) const +PointCloudBVH::allNearestNeighbors(std::size_t a_k) const { EBGEOMETRY_EXPECT(a_k >= 1); - EBGEOMETRY_EXPECT(a_eps >= T(0)); const std::size_t numPoints = m_positions.size(); - // Approximate-query pruning factor: a box is explored only if it could beat best/(1+eps), so any - // returned neighbor is within a factor (1+a_eps) of the true nearest. a_eps == 0 (the default) makes - // the factor 1 and the search exact. Distances are squared throughout, hence the squared reciprocal. - const T pruneScale = T(1) / ((T(1) + a_eps) * (T(1) + a_eps)); - std::vector result(numPoints * a_k); // Process points in leaf (build) order: the top-down build already grouped them spatially, so @@ -517,7 +500,7 @@ PointCloudBVH::allNearestNeighbors(std::size_t a_k, T a_eps) cons std::size_t found = 0; - this->query(m_positions[p], a_k, &result[p * a_k], found, p, m_leafOff[p], m_leafCnt[p], pruneScale); + this->query(m_positions[p], a_k, &result[p * a_k], found, p, m_leafOff[p], m_leafCnt[p]); } return result; diff --git a/Tests/TestPointCloudBVH.cpp b/Tests/TestPointCloudBVH.cpp index 7b1c43ac..21d4a6e9 100644 --- a/Tests/TestPointCloudBVH.cpp +++ b/Tests/TestPointCloudBVH.cpp @@ -146,31 +146,6 @@ TEMPLATE_TEST_CASE("PointCloudBVH queries match brute force", "[PointCloudBVH]", } } - SECTION("eps-approximate allNearestNeighbors stays within the (1+eps) distance guarantee") - { - const auto exact = bvh.allNearestNeighbors(1); - REQUIRE(exact.size() == n); - - for (const T eps : {T(0.25), T(1.0), T(4.0)}) { - const T ratioBound = (T(1) + eps) * (T(1) + eps); // guarantee is on distance, we compare squares - const auto approx = bvh.allNearestNeighbors(1, eps); - REQUIRE(approx.size() == n); - - for (std::size_t i = 0; i < n; i++) { - // The returned neighbor is a real point, so it is never closer than the true nearest, and no - // farther than a factor (1+eps): exactTrue <= approx <= (1+eps)^2 * exactTrue (squared). - CHECK(approx[i].distanceSquared >= exact[i].distanceSquared * (T(1) - looseMargin())); - CHECK(approx[i].distanceSquared <= exact[i].distanceSquared * ratioBound * (T(1) + looseMargin())); - } - } - - // eps == 0 must reproduce the exact search bit-for-bit (the default path). - const auto exact0 = bvh.allNearestNeighbors(1, T(0)); - for (std::size_t i = 0; i < n; i++) { - CHECK_THAT(exact0[i].distanceSquared, withinAbsT(exact[i].distanceSquared, tol)); - } - } - SECTION("accessors return the stored cloud data") { for (std::size_t i = 0; i < n; i += 313) { From 8fca5043fff6dbeadc2b8ca70d1d311d3f0cebe6 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 16 Jul 2026 07:42:22 +0000 Subject: [PATCH 08/17] Cap the default BVH branching ratio at 8 for float on AVX-512F (cache-line fit) DefaultBranchingRatio() returned 16 on AVX-512F (the SIMD-register-filling value). But the flat node stores K child offsets, so a K=16 float node is 96 B (24 B bounding volume + 8 B leaf fields + 16x4 B child offsets) and straddles two 64-byte cache lines, whereas K=8 is exactly 64 B -- one line. Nearest-neighbor traversal is memory-latency-bound, and the cache-line-sized node measured ~10-13% faster on PointCloudBVH all-nearest-neighbor (both uniform and sphere-surface distributions) than the wider K=16 fan-out, so cap float at 8 to match double. Every other (ISA, precision) default already fit one cache line, so only the AVX-512F/float case changes. The K=16/float SIMD path in pruneTraverse is kept for callers who request K=16 explicitly; it is simply no longer the default. Leaf SoA width (TriangleSoA/PointSoA DefaultWidth, still 16 for float) is unchanged. Update the DefaultBranchingRatio docs, the pruneTraverse configuration note, the Parser/PointCloudBVH @tparam comments, and the ImplemBVH/Parsers .rst pages (default-K table cell and rules of thumb) to reflect the cache-driven cap. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01JhdRJydutVJqQjg2Fuee9Z --- Docs/Sphinx/source/ImplemBVH.rst | 24 ++++++++++++-------- Docs/Sphinx/source/Parsers.rst | 7 +++--- Source/EBGeometry_BVH.hpp | 35 ++++++++++++++++++----------- Source/EBGeometry_BVHImplem.hpp | 7 +++++- Source/EBGeometry_Parser.hpp | 7 +++--- Source/EBGeometry_PointCloudBVH.hpp | 4 +++- 6 files changed, 54 insertions(+), 30 deletions(-) diff --git a/Docs/Sphinx/source/ImplemBVH.rst b/Docs/Sphinx/source/ImplemBVH.rst index 68c21b9c..6fbb5c8d 100644 --- a/Docs/Sphinx/source/ImplemBVH.rst +++ b/Docs/Sphinx/source/ImplemBVH.rst @@ -611,13 +611,16 @@ so its packed data is shared once per wrapper regardless of how that one object' ``Parser::readIntoTriangleBVH`` mirrors ``TriMeshSDF`` (``StoragePolicy`` as its 5th template parameter, defaulting to ``BVH::ValueStorage>``). See :ref:`Chap:Parsers`. -SIMD-optimal K and W by ISA +Default K and W by ISA ______________________________ -The helper ``BVH::DefaultBranchingRatio()`` returns the SIMD-optimal branching factor for -the current compilation target. ``EBGeometry::TriangleSoA::DefaultWidth()`` gives the -matching SoA width. Both are used as template defaults for ``TriMeshSDF`` and -``Parser::readIntoTriangleBVH``. +The helper ``BVH::DefaultBranchingRatio()`` returns the default branching factor for the current +compilation target, and ``EBGeometry::TriangleSoA::DefaultWidth()`` gives the matching SoA leaf +width. Both are used as template defaults for ``TriMeshSDF`` and ``Parser::readIntoTriangleBVH``. +The branching factor balances SIMD width against cache footprint: it is the SIMD-register-filling +value in every case *except* ``float`` on AVX-512F, where it is capped at 8 (not the register-filling +16) so the flat node stays within a single 64-byte cache line. Traversal is memory-latency-bound, and +that cache-line-sized node measured faster than the wider fan-out. The SoA leaf width is unaffected. .. list-table:: Default K and W by ISA and precision :widths: 25 25 25 25 @@ -629,7 +632,7 @@ matching SoA width. Both are used as template defaults for ``TriMeshSDF`` and - ``TriangleSoA::DefaultWidth()`` * - AVX-512F - ``float`` - - 16 + - 8 (capped for cache; see below) - 16 * - AVX-512F - ``double`` @@ -678,6 +681,9 @@ Rules of thumb: while the SAH/TopDown partitioner is still free to split down to smaller, tighter leaves wherever the geometry calls for it. A leaf smaller than ``W`` simply pads its SoA block's unused lanes. -* ``K = BVH::DefaultBranchingRatio()`` is a good default. With AVX-512F - available you can try ``K = 16`` (float) — the child-AABB test is evaluated in - a single SIMD batch, and the wider fan-out reduces tree depth. +* ``K = BVH::DefaultBranchingRatio()`` is a good default. On AVX-512F with + ``float`` it is 8 (not 16), so the flat node fits one 64-byte cache line; since + traversal is memory-latency-bound this beat the wider fan-out in benchmarks. You + can still request ``K = 16`` (float) explicitly — the child-AABB test then fits a + single 512-bit SIMD batch and the tree is shallower, at the cost of a node that + straddles two cache lines — but measure before assuming the wider fan-out wins. diff --git a/Docs/Sphinx/source/Parsers.rst b/Docs/Sphinx/source/Parsers.rst index 8455cc4f..ffb945e3 100644 --- a/Docs/Sphinx/source/Parsers.rst +++ b/Docs/Sphinx/source/Parsers.rst @@ -136,9 +136,10 @@ ________________________________ ``readIntoTriangleBVH(filename, maxLeafGroups, build)`` converts all DCEL polygons to triangles, packs them into SoA groups of ``W``, and builds a ``PackedBVH``, returning a ``shared_ptr>`` (or a vector thereof). SIMD -intrinsics evaluate up to ``W`` triangles per leaf visit. ``K`` and ``W`` default to the -SIMD-optimal values for ``T`` on the current ISA (``BVH::DefaultBranchingRatio()`` and -``TriangleSoA::DefaultWidth()``, see :ref:`Chap:MeshSDFClasses`); ``maxLeafGroups`` (default 4) +intrinsics evaluate up to ``W`` triangles per leaf visit. ``K`` and ``W`` default to +``BVH::DefaultBranchingRatio()`` and ``TriangleSoA::DefaultWidth()`` for ``T`` on the current +ISA (``K`` balances SIMD width against node cache footprint -- capped at 8 for ``float`` on AVX-512F +so the node fits one cache line; see :ref:`Chap:MeshSDFClasses`); ``maxLeafGroups`` (default 4) bounds the number of full ``W``-sized SoA groups per BVH leaf; ``StoragePolicy`` defaults to ``BVH::ValueStorage>``, matching ``TriMeshSDF``'s own default (see :ref:`Chap:MeshSDFClasses` for the rationale, and why ``readIntoPackedBVH``/``MeshSDF`` above has diff --git a/Source/EBGeometry_BVH.hpp b/Source/EBGeometry_BVH.hpp index ebb2b628..e31aa27b 100644 --- a/Source/EBGeometry_BVH.hpp +++ b/Source/EBGeometry_BVH.hpp @@ -71,20 +71,28 @@ struct ClusterSpec }; /** - * @brief Returns the SIMD-optimal BVH branching factor for type T on the current target ISA. - * @details Maps the floating-point type and the compile-time ISA to the K that fills one - * SIMD register exactly: + * @brief Returns the default BVH branching factor for type T on the current target ISA. + * @details Balances SIMD width against cache footprint. A wider K amortizes each traversal step over + * more children, but a node stores K child offsets, so a wider K spills the flat node across more + * cache lines. Traversal (e.g. PointCloudBVH nearest-neighbor) is memory-latency-bound, so the node + * fitting in a single cache line matters more than filling the SIMD register: * * | ISA | T=float | T=double | * |-----------|---------|----------| - * | AVX-512F | 16 | 8 | + * | AVX-512F | 8 | 8 | * | AVX | 8 | 4 | * | SSE4.1 | 4 | 4 | * | fallback | 4 | 4 | * + * The only case that is *not* register-filling is T=float on AVX-512F: the register-filling value + * there is 16, but a K=16 float node (24 B bounding volume + 8 B leaf fields + 16x4 B child offsets + * = 96 B) straddles two 64-byte cache lines, whereas K=8 (24 + 8 + 8x4 = 64 B) is exactly one. That + * cache-line-sized node measured faster in nearest-neighbor benchmarks than the wider SIMD fan-out, + * so float is capped at 8 to match double. Every other entry already fits one line. + * * Usage: `size_t K = BVH::DefaultBranchingRatio()` as a template-parameter default. * @tparam T Floating-point precision type (float or double). - * @return Optimal K for the current ISA and T. + * @return Default K for the current ISA and T. */ template [[nodiscard]] constexpr size_t @@ -92,12 +100,11 @@ DefaultBranchingRatio() noexcept { static_assert(std::is_floating_point_v, "BVH::DefaultBranchingRatio requires a floating-point T"); #if defined(__AVX512F__) - if constexpr (std::is_same_v) { - return 8; - } - else { - return 16; - } + // Both precisions use K=8. For double this is the SIMD-register-filling value. For float the + // register-filling value would be 16, but we deliberately cap at 8 so the flat node stays within a + // single 64-byte cache line instead of straddling two -- traversal is memory-latency-bound, and the + // cache-line-sized node beat the wider K=16 fan-out in nearest-neighbor benchmarks. + return 8; #elif defined(__AVX__) if constexpr (std::is_same_v) { return 4; @@ -694,8 +701,10 @@ SAHKWaySplit(PrimAndBVList& a_list, * into std::floor(K/2) and std::ceil(K/2) subsets — exact for power-of-two K; a reasonable * approximation for other values. * - * Recommended K values by ISA: - * - AVX-512F, float → K=16 (one @c _mm512_load_ps covers all K children) + * SIMD-register-filling K values by ISA (what one aligned load covers). These are the + * BVH::DefaultBranchingRatio() defaults too, except AVX-512F/float, whose default is capped at 8 + * (not 16) so the node fits one cache line -- see that helper's docs: + * - AVX-512F, float → K=16 (one @c _mm512_load_ps covers all K children; default is 8, not 16) * - AVX-512F, double → K=8 (one @c _mm512_load_pd covers all K children) * - AVX, float → K=8 (one @c _mm256_load_ps) * - AVX, double → K=4 (one @c _mm256_load_pd) diff --git a/Source/EBGeometry_BVHImplem.hpp b/Source/EBGeometry_BVHImplem.hpp index bc3849a0..96dcb80f 100644 --- a/Source/EBGeometry_BVHImplem.hpp +++ b/Source/EBGeometry_BVHImplem.hpp @@ -1120,7 +1120,12 @@ PackedBVH::pruneTraverse(const Vec3T& a_point, // for both (K=8, T=double) and (K=16, T=float). _mm512_load_pd / _mm512_load_ps // both require 64-byte alignment — the static_assert below catches any mismatch. // - // Recommended configurations on AVX-512 hardware: + // The K=16/float path below still fills a 512-bit register in one load and is kept for callers who + // opt into K=16 explicitly, but it is no longer the default: DefaultBranchingRatio() caps K + // at 8 on AVX-512F so the flat node fits one cache line (see its docs). With the default K=8, float + // BVHs take the AVX (_mm256, K=8/float) child-distance path further below, not this K=16 branch. + // + // AVX-512 register-filling configurations (K=16/float requires an explicit K): // float → K=16, W=16 (one _mm512_load_ps covers all children and one leaf group) // double → K=8, W=8 (one _mm512_load_pd covers all children; AVX-512F replaces // the 2×_mm256_load_pd emulation in the AVX fallback below) diff --git a/Source/EBGeometry_Parser.hpp b/Source/EBGeometry_Parser.hpp index 325fd398..c52bbdf9 100644 --- a/Source/EBGeometry_Parser.hpp +++ b/Source/EBGeometry_Parser.hpp @@ -238,9 +238,10 @@ readIntoPackedBVH(const std::vector& a_files, const BVH::Build a_bu * At query time the BVH uses SIMD intrinsics to evaluate W triangles per leaf visit. * @tparam T Floating-point precision for signed-distance evaluation. * @tparam Meta Per-face metadata type. - * @tparam K BVH branching factor. Defaults to BVH::DefaultBranchingRatio() — the SIMD-optimal value for - * T on the current ISA (K=16/float or K=8/double on AVX-512F; K=8/float or K=4/double - * on AVX; K=4 otherwise). Override only when benchmarking or using non-SIMD builds. + * @tparam K BVH branching factor. Defaults to BVH::DefaultBranchingRatio() — the SIMD/cache + * balanced value for T on the current ISA (K=8 for both float and double on AVX-512F, float capped at + * 8 rather than 16 so the node fits one cache line; K=8/float or K=4/double on AVX; K=4 otherwise). + * Override only when benchmarking or using non-SIMD builds. * @tparam W SIMD lane width: triangles per SoA group. Defaults to TriangleSoA::DefaultWidth() * (8/float or 4/double on AVX; 4 otherwise). * @tparam StoragePolicy PackedBVH primitive storage policy forwarded to TriMeshSDF (see diff --git a/Source/EBGeometry_PointCloudBVH.hpp b/Source/EBGeometry_PointCloudBVH.hpp index fd7c9b19..3376b160 100644 --- a/Source/EBGeometry_PointCloudBVH.hpp +++ b/Source/EBGeometry_PointCloudBVH.hpp @@ -50,7 +50,9 @@ namespace EBGeometry { * @tparam T Floating-point precision. * @tparam Meta User metadata type stored per point and returned via metadata(). Defaults to the * cloud index itself (std::size_t). - * @tparam K BVH branching factor. Defaults to the SIMD-optimal value for T. + * @tparam K BVH branching factor. Defaults to BVH::DefaultBranchingRatio(), which balances SIMD + * width against node cache footprint (float is capped at 8 on AVX-512F so the node fits + * one cache line -- see that helper). * @tparam W Points per SoA leaf lane group. Defaults to the SIMD-optimal width for T. */ template Date: Thu, 16 Jul 2026 08:02:07 +0000 Subject: [PATCH 09/17] Document the float K=8 branching-ratio finding in the benchmark README Add a "Finding" section to Benchmark/README.md recording why BVH::DefaultBranchingRatio() is capped at 8 on AVX-512: a K=8 float node is one 64-byte cache line (24 B BV + 8 B leaf fields + 8x4 B child offsets) vs two lines at K=16, and both traversals are memory-latency-bound, so the cache-line-sized node measured faster (~7-14% on NearestNeighbor, ~6-9% query plus ~8% build on MeshSDF) despite the narrower SIMD fan-out. Notes the mesh case gives up a real 512-bit-load advantage and still wins, and the corollary that global float alone does not help query (K=16 re-inflates the node). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01JhdRJydutVJqQjg2Fuee9Z --- Benchmark/README.md | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/Benchmark/README.md b/Benchmark/README.md index 43cb752d..42f4f63c 100644 --- a/Benchmark/README.md +++ b/Benchmark/README.md @@ -94,3 +94,32 @@ TriangleMeshDistance ~115 10.5 widely-used point of reference. - `TriMeshSDF` and TriangleMeshDistance compute the *signed* distance (their purpose); the comparison is on unsigned closest-surface distance, and the cross-check confirms all three agree. + +Finding: the `float` BVH branching ratio (cache lines over SIMD width) +---------------------------------------------------------------------- + +These benchmarks drove a change to EBGeometry's default BVH branching factor, +`BVH::DefaultBranchingRatio()`. On AVX-512 the register-filling factor for `float` is 16 (one +`_mm512_load_ps` covers all children), which used to be the default. But the flat BVH node stores `K` +child offsets, so a `K=16` float node is 96 bytes — 24 B bounding volume + 8 B leaf fields + 16×4 B +child offsets — and straddles **two** 64-byte cache lines, whereas `K=8` is exactly 64 B: **one** line. + +Both traversals here turn out to be **memory-latency-bound**, not compute-bound, so the +cache-line-sized node wins despite the narrower SIMD fan-out (measured, one machine, `K=8` vs `K=16` +at `float` precision): + +- **`NearestNeighbor` (point cloud):** ~7% faster query on the uniform cube, ~14% faster on the + sphere surface. +- **`MeshSDF` (triangle mesh):** ~6–9% faster query and ~8% faster build. + +The mesh case is the telling one: its external queries use the SIMD child-box traversal where `K=16` +genuinely fills a 512-bit register — a real compute advantage `K=8` gives up for a 256-bit load — and +`K=8` *still* won. The single-cache-line node beats the wider load. So `float` is now capped at 8 to +match `double`; every other (ISA, precision) default already fit one cache line. `K=16` remains +available by passing the template parameter explicitly for workloads that favor it. + +A corollary worth flagging: simply compiling everything in `float` does **not** speed the query up on +its own — `float`'s old default `K=16` re-inflated the node and offset the smaller scalars, so the two +effects roughly cancelled. Holding `K=8` is what captures the win (build time, separately, is ~35% +faster in `float` regardless). As always these are single-machine snapshots — the point is the +mechanism and that it is reproducible, not the exact percentages. From dcad7e15c64e97f7134699d960b050556c43ce40 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 16 Jul 2026 08:57:14 +0000 Subject: [PATCH 10/17] Add an experimental packet-traversal allNearestNeighbors prototype allNearestNeighborsPacket(packetSize) walks a packet of spatially-adjacent queries (consecutive in build order) through the tree together: each node is loaded once for the whole packet, a per-node uint64 active mask tracks which members still need the subtree, and each member seeds from its own leaf and excludes its own point. Results are bit-identical to serial allNearestNeighbors (new unit-test section checks packet sizes 1/3/8/64 against the serial graph for both precisions); the idea was to expose memory-level parallelism on the otherwise serial, latency-exposed dependent node-load chain. Measured outcome: it does not help this workload -- a consistent ~3-17% slower than serial at every packet size, and still slower even on a 4M-point cloud whose working set spills the 260 MB L3. The reason is that allNearestNeighbors already processes queries in spatial (m_order) order, so consecutive queries traverse an overlapping, slowly-sliding window of nodes/leaves that stays cache-resident regardless of total cloud size. The batch ordering already hides the latency packets aim to hide, so packets only add bookkeeping (mask iteration, leaf re-checks, fatter stack entries). Kept as an experimental method (clearly marked): the technique could still pay off for a batch of unordered external query points, where serial gets no cross-query cache reuse. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01JhdRJydutVJqQjg2Fuee9Z --- Source/EBGeometry_PointCloudBVH.hpp | 16 +++ Source/EBGeometry_PointCloudBVHImplem.hpp | 135 ++++++++++++++++++++++ Tests/TestPointCloudBVH.cpp | 17 +++ 3 files changed, 168 insertions(+) diff --git a/Source/EBGeometry_PointCloudBVH.hpp b/Source/EBGeometry_PointCloudBVH.hpp index 3376b160..a8fcb554 100644 --- a/Source/EBGeometry_PointCloudBVH.hpp +++ b/Source/EBGeometry_PointCloudBVH.hpp @@ -162,6 +162,22 @@ class PointCloudBVH [[nodiscard]] inline std::vector allNearestNeighbors(std::size_t a_k = 1) const; + /** + * @brief Experimental: packet variant of allNearestNeighbors(1) that walks a_packetSize spatially + * adjacent queries through the tree together. + * @details Same result as allNearestNeighbors(1) (each point's single nearest *other* point), but + * instead of one independent descent per point it processes a packet of consecutive (in build + * order, hence spatially near) points at once: every node is loaded once for the whole packet, and + * a per-node active bitmask tracks which packet members still need it. The single-query descent is + * a serial chain of dependent, latency-exposed node loads; a packet turns those into several + * independent chains in flight, exposing the memory-level parallelism that chain otherwise wastes. + * Each query still seeds from its own leaf and excludes its own point, so results are identical. + * @param[in] a_packetSize Queries processed together (clamped to [1, 64]). ~8 is a reasonable start. + * @return numPoints() Hits; result[i] is the nearest other point of cloud point i. + */ + [[nodiscard]] inline std::vector + allNearestNeighborsPacket(std::size_t a_packetSize = 8) const; + /** * @brief Brute-force closest point to an arbitrary query point (O(N) reference for closestPoint()). * @details Full linear scan. Same result contract as closestPoint(). diff --git a/Source/EBGeometry_PointCloudBVHImplem.hpp b/Source/EBGeometry_PointCloudBVHImplem.hpp index 25d6213b..70a32000 100644 --- a/Source/EBGeometry_PointCloudBVHImplem.hpp +++ b/Source/EBGeometry_PointCloudBVHImplem.hpp @@ -15,8 +15,11 @@ #include #include #include +#include +#include #include #include +#include // Our includes #include "EBGeometry_Macros.hpp" @@ -506,6 +509,138 @@ PointCloudBVH::allNearestNeighbors(std::size_t a_k) const return result; } +template +inline std::vector::Hit> +PointCloudBVH::allNearestNeighborsPacket(std::size_t a_packetSize) const +{ + const std::size_t numPoints = m_positions.size(); + + std::vector result(numPoints); + + if (this->m_linearNodes.empty() || numPoints == 0) { + return result; + } + + constexpr std::size_t maxPacket = 64; // one uint64_t active mask per stack entry + const std::size_t packet = std::min(std::max(a_packetSize, 1), maxPacket); + + // Per-query traversal state. The query point, its running best (squared distance + cloud index), + // and the cloud index to exclude (the query point itself). + struct QState + { + Vec3T point; + T bestDistanceSquared; + std::size_t bestIndex; + std::size_t exclude; + }; + + // Scan one leaf (a_cnt SoA groups from a_off) for a single query, updating its best. Identical in + // spirit to the single-query scanLeafBest, just parameterized on the QState. + const auto scanLeaf = [this](QState& a_q, std::uint32_t a_off, std::uint32_t a_cnt) noexcept { + for (std::uint32_t g = 0; g < a_cnt; g++) { + const PointGroup& group = this->m_primitives[a_off + g]; + const std::array distances = group.getDistances2(a_q.point); + + for (std::size_t lane = 0; lane < W; lane++) { + const std::size_t cloudIndex = group.getMetaData(lane); + + if (cloudIndex != a_q.exclude && distances[lane] < a_q.bestDistanceSquared) { + a_q.bestDistanceSquared = distances[lane]; + a_q.bestIndex = cloudIndex; + } + } + } + }; + + std::array q; + + struct StackEntry + { + std::uint32_t nodeIdx; + std::uint64_t active; // bit i set => packet member i still needs this subtree + }; + + constexpr int maxStack = 256; + + StackEntry stack[maxStack]; + + for (std::size_t base = 0; base < numPoints; base += packet) { + const std::size_t qn = std::min(packet, numPoints - base); + + // Initialize each packet member and seed its bound from its own leaf (the tight up-front bound the + // single-query path relies on; without it the shared unordered descent would explore far regions). + for (std::size_t i = 0; i < qn; i++) { + const std::uint32_t p = m_order[base + i]; + + q[i] = QState{m_positions[p], std::numeric_limits::max(), s_none, p}; + + scanLeaf(q[i], m_leafOff[p], m_leafCnt[p]); + } + + const std::uint64_t allActive = (qn == 64) ? ~std::uint64_t(0) : ((std::uint64_t(1) << qn) - 1); + + int stackTop = 0; + + stack[stackTop++] = StackEntry{0U, allActive}; + + while (stackTop > 0) { + const StackEntry entry = stack[--stackTop]; + const Node& node = this->m_linearNodes[entry.nodeIdx]; + + if (node.isLeaf()) { + const std::uint32_t primOffset = node.getPrimitivesOffset(); + const std::uint32_t primCount = node.getNumPrimitives(); + + // Scan this leaf once, reused across every still-relevant packet member; re-check each + // member's bound first, since it may have tightened below this node since it was pushed. + std::uint64_t m = entry.active; + while (m != 0) { + const int i = __builtin_ctzll(m); + m &= m - 1; + + if (node.getDistanceToBoundingVolume2(q[i].point) < q[i].bestDistanceSquared) { + scanLeaf(q[i], primOffset, primCount); + } + } + } + else { + const auto& childOffsets = node.getChildOffsets(); + + for (std::size_t k = 0; k < K; k++) { + const Node& child = this->m_linearNodes[childOffsets[k]]; + std::uint64_t childMask = 0; + + std::uint64_t m = entry.active; + while (m != 0) { + const int i = __builtin_ctzll(m); + m &= m - 1; + + if (child.getDistanceToBoundingVolume2(q[i].point) < q[i].bestDistanceSquared) { + childMask |= (std::uint64_t(1) << i); + } + } + + if (childMask != 0) { + EBGEOMETRY_EXPECT(stackTop < maxStack); + + stack[stackTop++] = StackEntry{childOffsets[k], childMask}; + } + } + } + } + + for (std::size_t i = 0; i < qn; i++) { + const std::uint32_t p = m_order[base + i]; + + if (q[i].bestIndex != s_none) { + result[p] = Hit{q[i].bestIndex, q[i].bestDistanceSquared}; + } + } + } + + return result; +} + template inline typename PointCloudBVH::Hit PointCloudBVH::bruteForceOne(const Vec3T& a_query, std::size_t a_exclude) const noexcept diff --git a/Tests/TestPointCloudBVH.cpp b/Tests/TestPointCloudBVH.cpp index 21d4a6e9..2b060fec 100644 --- a/Tests/TestPointCloudBVH.cpp +++ b/Tests/TestPointCloudBVH.cpp @@ -146,6 +146,23 @@ TEMPLATE_TEST_CASE("PointCloudBVH queries match brute force", "[PointCloudBVH]", } } + SECTION("allNearestNeighborsPacket matches serial allNearestNeighbors for every packet size") + { + const auto serial = bvh.allNearestNeighbors(1); + REQUIRE(serial.size() == n); + + // Packet sizes spanning single-query, sub-word, and full 64-bit mask; results must be identical + // regardless of how many queries share the walk. + for (const std::size_t packet : {std::size_t(1), std::size_t(3), std::size_t(8), std::size_t(64)}) { + const auto pk = bvh.allNearestNeighborsPacket(packet); + REQUIRE(pk.size() == n); + + for (std::size_t i = 0; i < n; i++) { + CHECK_THAT(pk[i].distanceSquared, withinAbsT(serial[i].distanceSquared, tol)); + } + } + } + SECTION("accessors return the stored cloud data") { for (std::size_t i = 0; i < n; i += 313) { From 844a849306ca360d594f680ac5cd29792e5cabb9 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 16 Jul 2026 09:22:32 +0000 Subject: [PATCH 11/17] Revert "Add an experimental packet-traversal allNearestNeighbors prototype" This reverts commit dcad7e15c64e97f7134699d960b050556c43ce40. --- Source/EBGeometry_PointCloudBVH.hpp | 16 --- Source/EBGeometry_PointCloudBVHImplem.hpp | 135 ---------------------- Tests/TestPointCloudBVH.cpp | 17 --- 3 files changed, 168 deletions(-) diff --git a/Source/EBGeometry_PointCloudBVH.hpp b/Source/EBGeometry_PointCloudBVH.hpp index a8fcb554..3376b160 100644 --- a/Source/EBGeometry_PointCloudBVH.hpp +++ b/Source/EBGeometry_PointCloudBVH.hpp @@ -162,22 +162,6 @@ class PointCloudBVH [[nodiscard]] inline std::vector allNearestNeighbors(std::size_t a_k = 1) const; - /** - * @brief Experimental: packet variant of allNearestNeighbors(1) that walks a_packetSize spatially - * adjacent queries through the tree together. - * @details Same result as allNearestNeighbors(1) (each point's single nearest *other* point), but - * instead of one independent descent per point it processes a packet of consecutive (in build - * order, hence spatially near) points at once: every node is loaded once for the whole packet, and - * a per-node active bitmask tracks which packet members still need it. The single-query descent is - * a serial chain of dependent, latency-exposed node loads; a packet turns those into several - * independent chains in flight, exposing the memory-level parallelism that chain otherwise wastes. - * Each query still seeds from its own leaf and excludes its own point, so results are identical. - * @param[in] a_packetSize Queries processed together (clamped to [1, 64]). ~8 is a reasonable start. - * @return numPoints() Hits; result[i] is the nearest other point of cloud point i. - */ - [[nodiscard]] inline std::vector - allNearestNeighborsPacket(std::size_t a_packetSize = 8) const; - /** * @brief Brute-force closest point to an arbitrary query point (O(N) reference for closestPoint()). * @details Full linear scan. Same result contract as closestPoint(). diff --git a/Source/EBGeometry_PointCloudBVHImplem.hpp b/Source/EBGeometry_PointCloudBVHImplem.hpp index 70a32000..25d6213b 100644 --- a/Source/EBGeometry_PointCloudBVHImplem.hpp +++ b/Source/EBGeometry_PointCloudBVHImplem.hpp @@ -15,11 +15,8 @@ #include #include #include -#include -#include #include #include -#include // Our includes #include "EBGeometry_Macros.hpp" @@ -509,138 +506,6 @@ PointCloudBVH::allNearestNeighbors(std::size_t a_k) const return result; } -template -inline std::vector::Hit> -PointCloudBVH::allNearestNeighborsPacket(std::size_t a_packetSize) const -{ - const std::size_t numPoints = m_positions.size(); - - std::vector result(numPoints); - - if (this->m_linearNodes.empty() || numPoints == 0) { - return result; - } - - constexpr std::size_t maxPacket = 64; // one uint64_t active mask per stack entry - const std::size_t packet = std::min(std::max(a_packetSize, 1), maxPacket); - - // Per-query traversal state. The query point, its running best (squared distance + cloud index), - // and the cloud index to exclude (the query point itself). - struct QState - { - Vec3T point; - T bestDistanceSquared; - std::size_t bestIndex; - std::size_t exclude; - }; - - // Scan one leaf (a_cnt SoA groups from a_off) for a single query, updating its best. Identical in - // spirit to the single-query scanLeafBest, just parameterized on the QState. - const auto scanLeaf = [this](QState& a_q, std::uint32_t a_off, std::uint32_t a_cnt) noexcept { - for (std::uint32_t g = 0; g < a_cnt; g++) { - const PointGroup& group = this->m_primitives[a_off + g]; - const std::array distances = group.getDistances2(a_q.point); - - for (std::size_t lane = 0; lane < W; lane++) { - const std::size_t cloudIndex = group.getMetaData(lane); - - if (cloudIndex != a_q.exclude && distances[lane] < a_q.bestDistanceSquared) { - a_q.bestDistanceSquared = distances[lane]; - a_q.bestIndex = cloudIndex; - } - } - } - }; - - std::array q; - - struct StackEntry - { - std::uint32_t nodeIdx; - std::uint64_t active; // bit i set => packet member i still needs this subtree - }; - - constexpr int maxStack = 256; - - StackEntry stack[maxStack]; - - for (std::size_t base = 0; base < numPoints; base += packet) { - const std::size_t qn = std::min(packet, numPoints - base); - - // Initialize each packet member and seed its bound from its own leaf (the tight up-front bound the - // single-query path relies on; without it the shared unordered descent would explore far regions). - for (std::size_t i = 0; i < qn; i++) { - const std::uint32_t p = m_order[base + i]; - - q[i] = QState{m_positions[p], std::numeric_limits::max(), s_none, p}; - - scanLeaf(q[i], m_leafOff[p], m_leafCnt[p]); - } - - const std::uint64_t allActive = (qn == 64) ? ~std::uint64_t(0) : ((std::uint64_t(1) << qn) - 1); - - int stackTop = 0; - - stack[stackTop++] = StackEntry{0U, allActive}; - - while (stackTop > 0) { - const StackEntry entry = stack[--stackTop]; - const Node& node = this->m_linearNodes[entry.nodeIdx]; - - if (node.isLeaf()) { - const std::uint32_t primOffset = node.getPrimitivesOffset(); - const std::uint32_t primCount = node.getNumPrimitives(); - - // Scan this leaf once, reused across every still-relevant packet member; re-check each - // member's bound first, since it may have tightened below this node since it was pushed. - std::uint64_t m = entry.active; - while (m != 0) { - const int i = __builtin_ctzll(m); - m &= m - 1; - - if (node.getDistanceToBoundingVolume2(q[i].point) < q[i].bestDistanceSquared) { - scanLeaf(q[i], primOffset, primCount); - } - } - } - else { - const auto& childOffsets = node.getChildOffsets(); - - for (std::size_t k = 0; k < K; k++) { - const Node& child = this->m_linearNodes[childOffsets[k]]; - std::uint64_t childMask = 0; - - std::uint64_t m = entry.active; - while (m != 0) { - const int i = __builtin_ctzll(m); - m &= m - 1; - - if (child.getDistanceToBoundingVolume2(q[i].point) < q[i].bestDistanceSquared) { - childMask |= (std::uint64_t(1) << i); - } - } - - if (childMask != 0) { - EBGEOMETRY_EXPECT(stackTop < maxStack); - - stack[stackTop++] = StackEntry{childOffsets[k], childMask}; - } - } - } - } - - for (std::size_t i = 0; i < qn; i++) { - const std::uint32_t p = m_order[base + i]; - - if (q[i].bestIndex != s_none) { - result[p] = Hit{q[i].bestIndex, q[i].bestDistanceSquared}; - } - } - } - - return result; -} - template inline typename PointCloudBVH::Hit PointCloudBVH::bruteForceOne(const Vec3T& a_query, std::size_t a_exclude) const noexcept diff --git a/Tests/TestPointCloudBVH.cpp b/Tests/TestPointCloudBVH.cpp index 2b060fec..21d4a6e9 100644 --- a/Tests/TestPointCloudBVH.cpp +++ b/Tests/TestPointCloudBVH.cpp @@ -146,23 +146,6 @@ TEMPLATE_TEST_CASE("PointCloudBVH queries match brute force", "[PointCloudBVH]", } } - SECTION("allNearestNeighborsPacket matches serial allNearestNeighbors for every packet size") - { - const auto serial = bvh.allNearestNeighbors(1); - REQUIRE(serial.size() == n); - - // Packet sizes spanning single-query, sub-word, and full 64-bit mask; results must be identical - // regardless of how many queries share the walk. - for (const std::size_t packet : {std::size_t(1), std::size_t(3), std::size_t(8), std::size_t(64)}) { - const auto pk = bvh.allNearestNeighborsPacket(packet); - REQUIRE(pk.size() == n); - - for (std::size_t i = 0; i < n; i++) { - CHECK_THAT(pk[i].distanceSquared, withinAbsT(serial[i].distanceSquared, tol)); - } - } - } - SECTION("accessors return the stored cloud data") { for (std::size_t i = 0; i < n; i += 313) { From 02e10773ecc83d6356de48cb3fb8ceea48ad6e4e Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 16 Jul 2026 18:08:32 +0000 Subject: [PATCH 12/17] Revert "Document the float K=8 branching-ratio finding in the benchmark README" This reverts commit 410f6893cdc56e6db9bf6a05aa40fe281df92361. --- Benchmark/README.md | 29 ----------------------------- 1 file changed, 29 deletions(-) diff --git a/Benchmark/README.md b/Benchmark/README.md index 42f4f63c..43cb752d 100644 --- a/Benchmark/README.md +++ b/Benchmark/README.md @@ -94,32 +94,3 @@ TriangleMeshDistance ~115 10.5 widely-used point of reference. - `TriMeshSDF` and TriangleMeshDistance compute the *signed* distance (their purpose); the comparison is on unsigned closest-surface distance, and the cross-check confirms all three agree. - -Finding: the `float` BVH branching ratio (cache lines over SIMD width) ----------------------------------------------------------------------- - -These benchmarks drove a change to EBGeometry's default BVH branching factor, -`BVH::DefaultBranchingRatio()`. On AVX-512 the register-filling factor for `float` is 16 (one -`_mm512_load_ps` covers all children), which used to be the default. But the flat BVH node stores `K` -child offsets, so a `K=16` float node is 96 bytes — 24 B bounding volume + 8 B leaf fields + 16×4 B -child offsets — and straddles **two** 64-byte cache lines, whereas `K=8` is exactly 64 B: **one** line. - -Both traversals here turn out to be **memory-latency-bound**, not compute-bound, so the -cache-line-sized node wins despite the narrower SIMD fan-out (measured, one machine, `K=8` vs `K=16` -at `float` precision): - -- **`NearestNeighbor` (point cloud):** ~7% faster query on the uniform cube, ~14% faster on the - sphere surface. -- **`MeshSDF` (triangle mesh):** ~6–9% faster query and ~8% faster build. - -The mesh case is the telling one: its external queries use the SIMD child-box traversal where `K=16` -genuinely fills a 512-bit register — a real compute advantage `K=8` gives up for a 256-bit load — and -`K=8` *still* won. The single-cache-line node beats the wider load. So `float` is now capped at 8 to -match `double`; every other (ISA, precision) default already fit one cache line. `K=16` remains -available by passing the template parameter explicitly for workloads that favor it. - -A corollary worth flagging: simply compiling everything in `float` does **not** speed the query up on -its own — `float`'s old default `K=16` re-inflated the node and offset the smaller scalars, so the two -effects roughly cancelled. Holding `K=8` is what captures the win (build time, separately, is ~35% -faster in `float` regardless). As always these are single-machine snapshots — the point is the -mechanism and that it is reproducible, not the exact percentages. From 823ad54719012220c7952e29f40ae7e1b04cc173 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 16 Jul 2026 18:08:32 +0000 Subject: [PATCH 13/17] Revert "Cap the default BVH branching ratio at 8 for float on AVX-512F (cache-line fit)" This reverts commit 8fca5043fff6dbeadc2b8ca70d1d311d3f0cebe6. --- Docs/Sphinx/source/ImplemBVH.rst | 24 ++++++++------------ Docs/Sphinx/source/Parsers.rst | 7 +++--- Source/EBGeometry_BVH.hpp | 35 +++++++++++------------------ Source/EBGeometry_BVHImplem.hpp | 7 +----- Source/EBGeometry_Parser.hpp | 7 +++--- Source/EBGeometry_PointCloudBVH.hpp | 4 +--- 6 files changed, 30 insertions(+), 54 deletions(-) diff --git a/Docs/Sphinx/source/ImplemBVH.rst b/Docs/Sphinx/source/ImplemBVH.rst index 6fbb5c8d..68c21b9c 100644 --- a/Docs/Sphinx/source/ImplemBVH.rst +++ b/Docs/Sphinx/source/ImplemBVH.rst @@ -611,16 +611,13 @@ so its packed data is shared once per wrapper regardless of how that one object' ``Parser::readIntoTriangleBVH`` mirrors ``TriMeshSDF`` (``StoragePolicy`` as its 5th template parameter, defaulting to ``BVH::ValueStorage>``). See :ref:`Chap:Parsers`. -Default K and W by ISA +SIMD-optimal K and W by ISA ______________________________ -The helper ``BVH::DefaultBranchingRatio()`` returns the default branching factor for the current -compilation target, and ``EBGeometry::TriangleSoA::DefaultWidth()`` gives the matching SoA leaf -width. Both are used as template defaults for ``TriMeshSDF`` and ``Parser::readIntoTriangleBVH``. -The branching factor balances SIMD width against cache footprint: it is the SIMD-register-filling -value in every case *except* ``float`` on AVX-512F, where it is capped at 8 (not the register-filling -16) so the flat node stays within a single 64-byte cache line. Traversal is memory-latency-bound, and -that cache-line-sized node measured faster than the wider fan-out. The SoA leaf width is unaffected. +The helper ``BVH::DefaultBranchingRatio()`` returns the SIMD-optimal branching factor for +the current compilation target. ``EBGeometry::TriangleSoA::DefaultWidth()`` gives the +matching SoA width. Both are used as template defaults for ``TriMeshSDF`` and +``Parser::readIntoTriangleBVH``. .. list-table:: Default K and W by ISA and precision :widths: 25 25 25 25 @@ -632,7 +629,7 @@ that cache-line-sized node measured faster than the wider fan-out. The SoA leaf - ``TriangleSoA::DefaultWidth()`` * - AVX-512F - ``float`` - - 8 (capped for cache; see below) + - 16 - 16 * - AVX-512F - ``double`` @@ -681,9 +678,6 @@ Rules of thumb: while the SAH/TopDown partitioner is still free to split down to smaller, tighter leaves wherever the geometry calls for it. A leaf smaller than ``W`` simply pads its SoA block's unused lanes. -* ``K = BVH::DefaultBranchingRatio()`` is a good default. On AVX-512F with - ``float`` it is 8 (not 16), so the flat node fits one 64-byte cache line; since - traversal is memory-latency-bound this beat the wider fan-out in benchmarks. You - can still request ``K = 16`` (float) explicitly — the child-AABB test then fits a - single 512-bit SIMD batch and the tree is shallower, at the cost of a node that - straddles two cache lines — but measure before assuming the wider fan-out wins. +* ``K = BVH::DefaultBranchingRatio()`` is a good default. With AVX-512F + available you can try ``K = 16`` (float) — the child-AABB test is evaluated in + a single SIMD batch, and the wider fan-out reduces tree depth. diff --git a/Docs/Sphinx/source/Parsers.rst b/Docs/Sphinx/source/Parsers.rst index ffb945e3..8455cc4f 100644 --- a/Docs/Sphinx/source/Parsers.rst +++ b/Docs/Sphinx/source/Parsers.rst @@ -136,10 +136,9 @@ ________________________________ ``readIntoTriangleBVH(filename, maxLeafGroups, build)`` converts all DCEL polygons to triangles, packs them into SoA groups of ``W``, and builds a ``PackedBVH``, returning a ``shared_ptr>`` (or a vector thereof). SIMD -intrinsics evaluate up to ``W`` triangles per leaf visit. ``K`` and ``W`` default to -``BVH::DefaultBranchingRatio()`` and ``TriangleSoA::DefaultWidth()`` for ``T`` on the current -ISA (``K`` balances SIMD width against node cache footprint -- capped at 8 for ``float`` on AVX-512F -so the node fits one cache line; see :ref:`Chap:MeshSDFClasses`); ``maxLeafGroups`` (default 4) +intrinsics evaluate up to ``W`` triangles per leaf visit. ``K`` and ``W`` default to the +SIMD-optimal values for ``T`` on the current ISA (``BVH::DefaultBranchingRatio()`` and +``TriangleSoA::DefaultWidth()``, see :ref:`Chap:MeshSDFClasses`); ``maxLeafGroups`` (default 4) bounds the number of full ``W``-sized SoA groups per BVH leaf; ``StoragePolicy`` defaults to ``BVH::ValueStorage>``, matching ``TriMeshSDF``'s own default (see :ref:`Chap:MeshSDFClasses` for the rationale, and why ``readIntoPackedBVH``/``MeshSDF`` above has diff --git a/Source/EBGeometry_BVH.hpp b/Source/EBGeometry_BVH.hpp index e31aa27b..ebb2b628 100644 --- a/Source/EBGeometry_BVH.hpp +++ b/Source/EBGeometry_BVH.hpp @@ -71,28 +71,20 @@ struct ClusterSpec }; /** - * @brief Returns the default BVH branching factor for type T on the current target ISA. - * @details Balances SIMD width against cache footprint. A wider K amortizes each traversal step over - * more children, but a node stores K child offsets, so a wider K spills the flat node across more - * cache lines. Traversal (e.g. PointCloudBVH nearest-neighbor) is memory-latency-bound, so the node - * fitting in a single cache line matters more than filling the SIMD register: + * @brief Returns the SIMD-optimal BVH branching factor for type T on the current target ISA. + * @details Maps the floating-point type and the compile-time ISA to the K that fills one + * SIMD register exactly: * * | ISA | T=float | T=double | * |-----------|---------|----------| - * | AVX-512F | 8 | 8 | + * | AVX-512F | 16 | 8 | * | AVX | 8 | 4 | * | SSE4.1 | 4 | 4 | * | fallback | 4 | 4 | * - * The only case that is *not* register-filling is T=float on AVX-512F: the register-filling value - * there is 16, but a K=16 float node (24 B bounding volume + 8 B leaf fields + 16x4 B child offsets - * = 96 B) straddles two 64-byte cache lines, whereas K=8 (24 + 8 + 8x4 = 64 B) is exactly one. That - * cache-line-sized node measured faster in nearest-neighbor benchmarks than the wider SIMD fan-out, - * so float is capped at 8 to match double. Every other entry already fits one line. - * * Usage: `size_t K = BVH::DefaultBranchingRatio()` as a template-parameter default. * @tparam T Floating-point precision type (float or double). - * @return Default K for the current ISA and T. + * @return Optimal K for the current ISA and T. */ template [[nodiscard]] constexpr size_t @@ -100,11 +92,12 @@ DefaultBranchingRatio() noexcept { static_assert(std::is_floating_point_v, "BVH::DefaultBranchingRatio requires a floating-point T"); #if defined(__AVX512F__) - // Both precisions use K=8. For double this is the SIMD-register-filling value. For float the - // register-filling value would be 16, but we deliberately cap at 8 so the flat node stays within a - // single 64-byte cache line instead of straddling two -- traversal is memory-latency-bound, and the - // cache-line-sized node beat the wider K=16 fan-out in nearest-neighbor benchmarks. - return 8; + if constexpr (std::is_same_v) { + return 8; + } + else { + return 16; + } #elif defined(__AVX__) if constexpr (std::is_same_v) { return 4; @@ -701,10 +694,8 @@ SAHKWaySplit(PrimAndBVList& a_list, * into std::floor(K/2) and std::ceil(K/2) subsets — exact for power-of-two K; a reasonable * approximation for other values. * - * SIMD-register-filling K values by ISA (what one aligned load covers). These are the - * BVH::DefaultBranchingRatio() defaults too, except AVX-512F/float, whose default is capped at 8 - * (not 16) so the node fits one cache line -- see that helper's docs: - * - AVX-512F, float → K=16 (one @c _mm512_load_ps covers all K children; default is 8, not 16) + * Recommended K values by ISA: + * - AVX-512F, float → K=16 (one @c _mm512_load_ps covers all K children) * - AVX-512F, double → K=8 (one @c _mm512_load_pd covers all K children) * - AVX, float → K=8 (one @c _mm256_load_ps) * - AVX, double → K=4 (one @c _mm256_load_pd) diff --git a/Source/EBGeometry_BVHImplem.hpp b/Source/EBGeometry_BVHImplem.hpp index 96dcb80f..bc3849a0 100644 --- a/Source/EBGeometry_BVHImplem.hpp +++ b/Source/EBGeometry_BVHImplem.hpp @@ -1120,12 +1120,7 @@ PackedBVH::pruneTraverse(const Vec3T& a_point, // for both (K=8, T=double) and (K=16, T=float). _mm512_load_pd / _mm512_load_ps // both require 64-byte alignment — the static_assert below catches any mismatch. // - // The K=16/float path below still fills a 512-bit register in one load and is kept for callers who - // opt into K=16 explicitly, but it is no longer the default: DefaultBranchingRatio() caps K - // at 8 on AVX-512F so the flat node fits one cache line (see its docs). With the default K=8, float - // BVHs take the AVX (_mm256, K=8/float) child-distance path further below, not this K=16 branch. - // - // AVX-512 register-filling configurations (K=16/float requires an explicit K): + // Recommended configurations on AVX-512 hardware: // float → K=16, W=16 (one _mm512_load_ps covers all children and one leaf group) // double → K=8, W=8 (one _mm512_load_pd covers all children; AVX-512F replaces // the 2×_mm256_load_pd emulation in the AVX fallback below) diff --git a/Source/EBGeometry_Parser.hpp b/Source/EBGeometry_Parser.hpp index c52bbdf9..325fd398 100644 --- a/Source/EBGeometry_Parser.hpp +++ b/Source/EBGeometry_Parser.hpp @@ -238,10 +238,9 @@ readIntoPackedBVH(const std::vector& a_files, const BVH::Build a_bu * At query time the BVH uses SIMD intrinsics to evaluate W triangles per leaf visit. * @tparam T Floating-point precision for signed-distance evaluation. * @tparam Meta Per-face metadata type. - * @tparam K BVH branching factor. Defaults to BVH::DefaultBranchingRatio() — the SIMD/cache - * balanced value for T on the current ISA (K=8 for both float and double on AVX-512F, float capped at - * 8 rather than 16 so the node fits one cache line; K=8/float or K=4/double on AVX; K=4 otherwise). - * Override only when benchmarking or using non-SIMD builds. + * @tparam K BVH branching factor. Defaults to BVH::DefaultBranchingRatio() — the SIMD-optimal value for + * T on the current ISA (K=16/float or K=8/double on AVX-512F; K=8/float or K=4/double + * on AVX; K=4 otherwise). Override only when benchmarking or using non-SIMD builds. * @tparam W SIMD lane width: triangles per SoA group. Defaults to TriangleSoA::DefaultWidth() * (8/float or 4/double on AVX; 4 otherwise). * @tparam StoragePolicy PackedBVH primitive storage policy forwarded to TriMeshSDF (see diff --git a/Source/EBGeometry_PointCloudBVH.hpp b/Source/EBGeometry_PointCloudBVH.hpp index 3376b160..fd7c9b19 100644 --- a/Source/EBGeometry_PointCloudBVH.hpp +++ b/Source/EBGeometry_PointCloudBVH.hpp @@ -50,9 +50,7 @@ namespace EBGeometry { * @tparam T Floating-point precision. * @tparam Meta User metadata type stored per point and returned via metadata(). Defaults to the * cloud index itself (std::size_t). - * @tparam K BVH branching factor. Defaults to BVH::DefaultBranchingRatio(), which balances SIMD - * width against node cache footprint (float is capped at 8 on AVX-512F so the node fits - * one cache line -- see that helper). + * @tparam K BVH branching factor. Defaults to the SIMD-optimal value for T. * @tparam W Points per SoA leaf lane group. Defaults to the SIMD-optimal width for T. */ template Date: Thu, 16 Jul 2026 20:02:59 +0000 Subject: [PATCH 14/17] Add kd3 to the NearestNeighbor benchmark comparison Add https://github.com/KaruroChori/kd3 (a header-only SoA/SIMD KD-tree for low-dimensional point kNN) as a Benchmark/kd3 submodule and a fifth method in the all-nearest-neighbor benchmark, alongside PointCloudBVH, PointCloudHashGrid, picoflann, and nanoflann. Each point's nearest other point is found via kd3's query_knn(k=2) in Hilbert order (skipping the query point itself), cross-checked against the same brute-force baseline as the others. kd3 requires C++23 (std::expected/std::span), so the whole NearestNeighbor benchmark is now built with -std=c++23 (EBGeometry/nanoflann/picoflann all compile unchanged under it). For a like-for-like comparison kd3 is run in double and single-threaded (compiled without -fopenmp); its distance_t is pinned to the benchmark's T rather than its float default. Its headline float/SIMD/OpenMP fast path is therefore not exercised here -- noted in the code and README so the numbers aren't misread. In this restricted mode kd3 posts the fastest build and a competitive query. Wire kd3 into the compile-only Benchmark-Compile CI job's submodule list, and document it in Benchmark/README.md and Docs/Sphinx/source/Benchmark.rst (library list, fetch command, results table, and the parity caveats). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01JhdRJydutVJqQjg2Fuee9Z --- .github/workflows/CI.yml | 2 +- .gitmodules | 4 ++ Benchmark/NearestNeighbor/GNUmakefile | 13 +++-- Benchmark/NearestNeighbor/main.cpp | 73 +++++++++++++++++++++++++-- Benchmark/README.md | 22 +++++--- Benchmark/kd3 | 1 + Docs/Sphinx/source/Benchmark.rst | 10 ++-- 7 files changed, 103 insertions(+), 22 deletions(-) create mode 160000 Benchmark/kd3 diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 41c5c3f6..7891cda9 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -445,7 +445,7 @@ jobs: matrix: include: - directory: Benchmark/NearestNeighbor - submodules: Benchmark/picoflann Benchmark/nanoflann + submodules: Benchmark/picoflann Benchmark/nanoflann Benchmark/kd3 - directory: Benchmark/MeshSDF submodules: Benchmark/fcpw Benchmark/TriangleMeshDistance steps: diff --git a/.gitmodules b/.gitmodules index 8a654ef1..486b857a 100644 --- a/.gitmodules +++ b/.gitmodules @@ -17,3 +17,7 @@ [submodule "Benchmark/TriangleMeshDistance"] path = Benchmark/TriangleMeshDistance url = https://github.com/InteractiveComputerGraphics/TriangleMeshDistance +[submodule "Benchmark/kd3"] + path = Benchmark/kd3 + url = https://github.com/KaruroChori/kd3 + shallow = true diff --git a/Benchmark/NearestNeighbor/GNUmakefile b/Benchmark/NearestNeighbor/GNUmakefile index 0007fbe9..54ecf841 100644 --- a/Benchmark/NearestNeighbor/GNUmakefile +++ b/Benchmark/NearestNeighbor/GNUmakefile @@ -1,15 +1,18 @@ -# Benchmark: PointCloudBVH vs picoflann vs nanoflann (all-nearest-neighbor). -# Requires the Benchmark/{picoflann,nanoflann} submodules: -# git submodule update --init Benchmark/picoflann Benchmark/nanoflann +# Benchmark: PointCloudBVH vs picoflann vs nanoflann vs kd3 (all-nearest-neighbor). +# Requires the Benchmark/{picoflann,nanoflann,kd3} submodules: +# git submodule update --init Benchmark/picoflann Benchmark/nanoflann Benchmark/kd3 +# kd3 needs C++23 (std::expected/std::span), so the whole benchmark is built with -std=c++23; it is +# compiled WITHOUT -fopenmp so kd3 (like the others) runs single-threaded, for a fair comparison. # EBGEOMETRY_HOME defaults to the repo root (two levels up). EBGEOMETRY_HOME ?= ../.. CXX ?= g++ -CXXFLAGS ?= -std=c++17 -O3 -march=native +CXXFLAGS ?= -std=c++23 -O3 -march=native INCLUDES := -I$(EBGEOMETRY_HOME) \ -I$(EBGEOMETRY_HOME)/Benchmark/picoflann \ - -I$(EBGEOMETRY_HOME)/Benchmark/nanoflann/include + -I$(EBGEOMETRY_HOME)/Benchmark/nanoflann/include \ + -I$(EBGEOMETRY_HOME)/Benchmark/kd3/include TARGET := NearestNeighbor.ex diff --git a/Benchmark/NearestNeighbor/main.cpp b/Benchmark/NearestNeighbor/main.cpp index 69bc1153..8d9eadfa 100644 --- a/Benchmark/NearestNeighbor/main.cpp +++ b/Benchmark/NearestNeighbor/main.cpp @@ -2,17 +2,23 @@ // // SPDX-License-Identifier: GPL-3.0-or-later -// Benchmark: EBGeometry PointCloudBVH and PointCloudHashGrid vs picoflann vs nanoflann, -// all-nearest-neighbor. 500k 3D points (double). Every point's nearest OTHER point. All four -// verified against a brute-force sample. Both KD-trees used vanilla; distances are squared throughout. +// Benchmark: EBGeometry PointCloudBVH and PointCloudHashGrid vs picoflann vs nanoflann vs kd3, +// all-nearest-neighbor. 500k 3D points (double). Every point's nearest OTHER point. All five +// verified against a brute-force sample. The KD-trees are used vanilla; distances are squared throughout. // // Two point distributions are benchmarked in turn, since spatial data structures behave very // differently depending on how the points fill space: // 1. Uniform in the unit cube -- points fill a 3D volume evenly (the easy, balanced case). // 2. On the unit-sphere surface -- points lie on a 2D manifold: locally dense, globally hollow, // a harder case for uniform grids (many empty cells inside). +// +// kd3 (https://github.com/KaruroChori/kd3) requires C++23 (std::expected/std::span), so this whole +// benchmark is built with -std=c++23. For a like-for-like comparison it is run here in double and +// single-threaded (compiled without -fopenmp); kd3 is primarily float/SIMD-tuned and can parallelize +// its build with OpenMP, neither of which is exercised here. #include +#include #include #include #include @@ -25,12 +31,17 @@ #include +#include + #include "nanoflann.hpp" #include "picoflann.h" using T = double; using Vec3 = EBGeometry::Vec3T; +// kd3 tree, in double precision (its distance_t defaults to float; pin it to T for a fair comparison). +using Kd3 = kd3::KdTree>; + constexpr std::size_t numPoints = 500000; constexpr std::size_t sampleSize = 500; constexpr std::uint64_t pointSeed = 123456789ULL; @@ -264,7 +275,59 @@ runCase(const std::string& a_label, const std::vector& a_positions) row("nanoflann", buildMs, queryUs, bad); } - std::cout << " flann query loops iterate in Hilbert order (warm cache, like EBGeometry's leaf order).\n"; + // ── kd3 (per-point query_knn, k=2) ── + { + // kd3 builds from its own FatPoint {coords, payload} array (payload = cloud index), and sorts it + // in place. Constructing that array is kd3's required input format, so it is folded into the timed + // build (mirroring how fcpw's soup construction is timed as part of its build in MeshSDF). + std::vector fat(n); + timer.start(); + for (std::size_t i = 0; i < n; i++) { + fat[i] = Kd3::FatPoint{{a_positions[i][0], a_positions[i][1], a_positions[i][2]}, static_cast(i)}; + } + + auto treeExpected = Kd3::build(fat); + timer.stop(); + const double buildMs = 1.0e3 * timer.seconds(); + + if (!treeExpected) { + row("kd3", buildMs, std::numeric_limits::infinity(), sampleSize); // build failed -> flag it + } + else { + const Kd3& tree = *treeExpected; + + auto nnOther = [&](std::size_t i) -> T { + const Kd3::point_t q = {a_positions[i][0], a_positions[i][1], a_positions[i][2]}; + std::array buf{}; + const auto res = tree.query_knn(q, buf); + if (res) { + for (const auto& kr : *res) { + if (kr.payload_id != static_cast(i)) { // skip the query point itself (dist 0) + return kr.dist_sq; + } + } + } + return std::numeric_limits::max(); + }; + + volatile T sink = T(0); + timer.start(); + for (const std::uint32_t p : order) { + sink += nnOther(p); + } + timer.stop(); + (void)sink; + const double queryUs = 1.0e6 * timer.seconds() / double(n); + + std::size_t bad = 0; + for (std::size_t s = 0; s < sampleSize; s++) { + bad += !ok(nnOther(s * stride), s); + } + row("kd3", buildMs, queryUs, bad); + } + } + + std::cout << " flann/kd3 query loops iterate in Hilbert order (warm cache, like EBGeometry's leaf order).\n"; std::cout << " One-time Hilbert sort the flann libs need for that order: " << std::setprecision(3) << sortUsPerPt << " us/pt\n (add to their query if counted; EBGeometry reuses its build order for free).\n\n"; } @@ -274,7 +337,7 @@ runCase(const std::string& a_label, const std::vector& a_positions) int main() { - std::cout << "All-nearest-neighbor: PointCloudBVH & PointCloudHashGrid vs picoflann vs nanoflann\n\n"; + std::cout << "All-nearest-neighbor: PointCloudBVH & PointCloudHashGrid vs picoflann vs nanoflann vs kd3\n\n"; runCase("Uniform in the unit cube", EBGeometry::Random::samplePoints(numPoints, pointSeed)); runCase("On the unit-sphere surface", samplePointsOnSphere(numPoints, pointSeed)); diff --git a/Benchmark/README.md b/Benchmark/README.md index 43cb752d..12f01b45 100644 --- a/Benchmark/README.md +++ b/Benchmark/README.md @@ -14,13 +14,14 @@ does not *run* or time them (see issue #109). They are pinned as git submodules * [nanoflann](https://github.com/jlblancoc/nanoflann) — header-only KD-tree (point kNN) * [picoflann](https://github.com/rmsalinas/picoflann) — tiny header-only KD-tree (point kNN) +* [kd3](https://github.com/KaruroChori/kd3) — header-only SoA/SIMD KD-tree (point kNN), **requires C++23** * [fcpw](https://github.com/rohan-sawhney/fcpw) — closest-point / SDF on triangle meshes * [TriangleMeshDistance](https://github.com/InteractiveComputerGraphics/TriangleMeshDistance) — header-only signed distance to triangle meshes Fetch them (and the top-level mesh submodule) with: ```bash -git submodule update --init Benchmark/nanoflann Benchmark/picoflann Benchmark/fcpw \ +git submodule update --init Benchmark/nanoflann Benchmark/picoflann Benchmark/kd3 Benchmark/fcpw \ Benchmark/TriangleMeshDistance common-3d-test-models git -C Benchmark/fcpw submodule update --init deps/eigen # fcpw's Eigen (skip the GPU slang-rhi dep) ``` @@ -36,7 +37,7 @@ that transfers across machines even when the timings do not. ---------------------------------------------------------- For every point in a 500,000-point cloud (double precision), find its nearest *other* point. Compares -EBGeometry's `PointCloudBVH` and `PointCloudHashGrid` against picoflann and nanoflann. The same +EBGeometry's `PointCloudBVH` and `PointCloudHashGrid` against picoflann, nanoflann, and kd3. The same comparison is run over two point distributions, since spatial structures behave very differently depending on how the points fill space: @@ -49,10 +50,11 @@ Representative result (one machine, illustrative — see the note on machine dep ``` uniform cube sphere surface Method Build(ms) Query(us/pt) Build(ms) Query(us/pt) -PointCloudBVH ~165 0.90 ~155 0.66 -PointCloudHashGrid ~19 1.85 ~17 2.04 -picoflann ~132 0.87 ~130 0.55 -nanoflann ~262 0.44 ~273 0.33 +PointCloudBVH ~118 0.62 ~105 0.48 +PointCloudHashGrid ~12 1.50 ~11 1.60 +picoflann ~104 0.75 ~104 0.48 +nanoflann ~225 0.41 ~230 0.36 +kd3 ~83 0.57 ~86 0.50 ``` - The KD-tree queries are iterated in **Hilbert order** so their node cache is as warm as @@ -61,13 +63,19 @@ nanoflann ~262 0.44 ~273 0.33 the EBGeometry structures reuse the ordering their build already produced. - **`PointCloudHashGrid` trades query speed for build speed**: an O(N) uniform grid builds ~8x faster than the BVH but scans neighbor cells per query, so it queries ~2x slower. It is the weakest on the - sphere surface (query ~2.0 us/pt) — the hollow distribution leaves its grid mostly empty while the + sphere surface (query ~1.6 us/pt) — the hollow distribution leaves its grid mostly empty while the occupied surface cells are denser than the ~1-point-per-cell target. - The tree/BVH methods, by contrast, get *faster* on the sphere surface than in the cube (the local neighborhood is effectively lower-dimensional, so pruning is tighter). nanoflann has the fastest raw per-query traversal throughout; `PointCloudBVH` is competitive end-to-end because it gets its query order for free and builds faster. Which structure to pick depends on the build/query balance and the point distribution — that is the point of running both cases. +- **kd3 is run here in a deliberately restricted mode for parity**: double precision and + single-threaded (compiled without `-fopenmp`). In that mode it posts the fastest *build* and a + competitive query. Its headline ">2x faster than nanoflann" claim, though, is for `float` with its + SoA/SIMD fast path and an OpenMP-parallel build — none of which this same-precision, same-thread + comparison exercises, so treat kd3's numbers here as a lower bound on what it can do. (kd3 requires + C++23, so this benchmark is built with `-std=c++23`.) `MeshSDF/` — closest-point on a triangle mesh --------------------------------------------- diff --git a/Benchmark/kd3 b/Benchmark/kd3 new file mode 160000 index 00000000..e389a4e5 --- /dev/null +++ b/Benchmark/kd3 @@ -0,0 +1 @@ +Subproject commit e389a4e51838f2dc3a367092013b11e91a3f6910 diff --git a/Docs/Sphinx/source/Benchmark.rst b/Docs/Sphinx/source/Benchmark.rst index 699b1b1e..2a79cb5b 100644 --- a/Docs/Sphinx/source/Benchmark.rst +++ b/Docs/Sphinx/source/Benchmark.rst @@ -20,16 +20,18 @@ submodule) with: .. code-block:: bash - git submodule update --init Benchmark/nanoflann Benchmark/picoflann Benchmark/fcpw \ + git submodule update --init Benchmark/nanoflann Benchmark/picoflann Benchmark/kd3 Benchmark/fcpw \ Benchmark/TriangleMeshDistance common-3d-test-models git -C Benchmark/fcpw submodule update --init deps/eigen # fcpw's Eigen (skip the GPU dep) Each benchmark ships a ``GNUmakefile`` (``make && ./.ex``). See each folder's ``README.md`` for the full detail and representative numbers. -* :file:`Benchmark/NearestNeighbor` -- all-nearest-neighbor over a point cloud: ``PointCloudBVH`` vs - `nanoflann `_ vs - `picoflann `_. +* :file:`Benchmark/NearestNeighbor` -- all-nearest-neighbor over a point cloud: ``PointCloudBVH`` and + ``PointCloudHashGrid`` vs `nanoflann `_ vs + `picoflann `_ vs + `kd3 `_ (a SoA/SIMD KD-tree; needs C++23, run here in double and + single-threaded for parity). * :file:`Benchmark/MeshSDF` -- closest-point on a triangle mesh: ``TriMeshSDF`` vs `fcpw `_ (built with its Enoki CPU vectorization) vs `TriangleMeshDistance `_. From 6546de5f795c66c06b7b9b4e51fddfd1c763f4cb Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 16 Jul 2026 20:05:55 +0000 Subject: [PATCH 15/17] Correct the kd3 caveat: its 2x query claim is single-threaded, not OpenMP The prior note wrongly implied kd3's "~2.2x vs nanoflann" query claim relied on an OpenMP-parallel build. It does not: in kd3's own benchmark both the kd3 and nanoflann query loops are plain serial loops (query_knn_inline is a noexcept inline call; the printed "parallelism: 32" is cosmetic), so the speedup is kd3's SoA/SIMD per-query kernel in float, single-threaded. OpenMP only parallelizes the build (a separate, disclosed claim). The real reason kd3 doesn't reach 2x in this benchmark is that it is run in double, halving its SIMD width -- not threading. Compiling without -fopenmp keeps the build comparison single-threaded and thus fairer than kd3's own OpenMP-parallel build numbers. Fix the wording in main.cpp, Benchmark/README.md, and Benchmark.rst accordingly. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01JhdRJydutVJqQjg2Fuee9Z --- Benchmark/NearestNeighbor/main.cpp | 10 +++++++--- Benchmark/README.md | 15 +++++++++------ Docs/Sphinx/source/Benchmark.rst | 4 ++-- 3 files changed, 18 insertions(+), 11 deletions(-) diff --git a/Benchmark/NearestNeighbor/main.cpp b/Benchmark/NearestNeighbor/main.cpp index 8d9eadfa..319052fe 100644 --- a/Benchmark/NearestNeighbor/main.cpp +++ b/Benchmark/NearestNeighbor/main.cpp @@ -13,9 +13,13 @@ // a harder case for uniform grids (many empty cells inside). // // kd3 (https://github.com/KaruroChori/kd3) requires C++23 (std::expected/std::span), so this whole -// benchmark is built with -std=c++23. For a like-for-like comparison it is run here in double and -// single-threaded (compiled without -fopenmp); kd3 is primarily float/SIMD-tuned and can parallelize -// its build with OpenMP, neither of which is exercised here. +// benchmark is built with -std=c++23. kd3's headline "~2.2x query throughput vs nanoflann" is a +// single-threaded FLOAT result (in kd3's own benchmark both query loops are plain serial loops; the +// speedup is its SoA/SIMD per-query kernel, not multithreading). Here kd3 is run in double, which +// halves that SIMD width -- the main reason it does not reach 2x in this table. It is also compiled +// without -fopenmp, so its build is single-threaded like everything else here; kd3's docs quote a +// much faster build, but that number is OpenMP-parallel (its author notes it is still faster serially). +// So these kd3 numbers are a fair single-threaded comparison, understated only by the double precision. #include #include diff --git a/Benchmark/README.md b/Benchmark/README.md index 12f01b45..923197b9 100644 --- a/Benchmark/README.md +++ b/Benchmark/README.md @@ -70,12 +70,15 @@ kd3 ~83 0.57 ~86 0.50 raw per-query traversal throughout; `PointCloudBVH` is competitive end-to-end because it gets its query order for free and builds faster. Which structure to pick depends on the build/query balance and the point distribution — that is the point of running both cases. -- **kd3 is run here in a deliberately restricted mode for parity**: double precision and - single-threaded (compiled without `-fopenmp`). In that mode it posts the fastest *build* and a - competitive query. Its headline ">2x faster than nanoflann" claim, though, is for `float` with its - SoA/SIMD fast path and an OpenMP-parallel build — none of which this same-precision, same-thread - comparison exercises, so treat kd3's numbers here as a lower bound on what it can do. (kd3 requires - C++23, so this benchmark is built with `-std=c++23`.) +- **kd3's headline "~2.2x query throughput vs nanoflann" is a single-threaded `float` result** — in + kd3's own benchmark both query loops are plain serial loops, so the speedup is its SoA/SIMD + per-query kernel, *not* multithreading. Here kd3 is run in **double**, which halves that SIMD width; + that is the main reason it doesn't reach 2x in this table, and it's why these numbers are a lower + bound on what kd3 can do. It is compiled **without `-fopenmp`**, so its build is single-threaded + like the others — a *fairer* build comparison than kd3's own docs, whose fast build number is + OpenMP-parallel (its author notes it is still faster serially). In this fair single-threaded mode + kd3 posts the fastest build and a competitive query. (kd3 needs C++23, so the benchmark is built + with `-std=c++23`.) `MeshSDF/` — closest-point on a triangle mesh --------------------------------------------- diff --git a/Docs/Sphinx/source/Benchmark.rst b/Docs/Sphinx/source/Benchmark.rst index 2a79cb5b..ba38bdb2 100644 --- a/Docs/Sphinx/source/Benchmark.rst +++ b/Docs/Sphinx/source/Benchmark.rst @@ -30,8 +30,8 @@ the full detail and representative numbers. * :file:`Benchmark/NearestNeighbor` -- all-nearest-neighbor over a point cloud: ``PointCloudBVH`` and ``PointCloudHashGrid`` vs `nanoflann `_ vs `picoflann `_ vs - `kd3 `_ (a SoA/SIMD KD-tree; needs C++23, run here in double and - single-threaded for parity). + `kd3 `_ (a SoA/SIMD KD-tree; needs C++23, run here in double -- + halving its float SIMD width -- for a same-precision comparison). * :file:`Benchmark/MeshSDF` -- closest-point on a triangle mesh: ``TriMeshSDF`` vs `fcpw `_ (built with its Enoki CPU vectorization) vs `TriangleMeshDistance `_. From 470dc18d57c436ed598c3bfd9fa3fbfd1c7bfabb Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 16 Jul 2026 20:10:58 +0000 Subject: [PATCH 16/17] Add a kd3 (float) row showing kd3's native SoA/SIMD best case Alongside the same-precision kd3 (double) row, run kd3 in its native float precision (kd3::limits) as a best-case reference -- its intended SoA/SIMD fast path -- and label the two rows kd3 (double) / kd3 (float). This is a different-precision reference, not apples-to-apples with the double field, mirroring how TriangleMeshDistance is reported (double) against float methods in the MeshSDF benchmark. The float row uses a looser cross-check tolerance since its squared distance is computed in float against the double brute-force truth. On this machine kd3 (float) runs ~7-11% faster than kd3 (double); the gain is modest here because AVX-512 clocks down, damping the float SIMD-width advantage. Update Benchmark/README.md (table + note) and Benchmark.rst accordingly. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01JhdRJydutVJqQjg2Fuee9Z --- Benchmark/NearestNeighbor/main.cpp | 68 ++++++++++++++++++++++++++++-- Benchmark/README.md | 9 +++- Docs/Sphinx/source/Benchmark.rst | 4 +- 3 files changed, 75 insertions(+), 6 deletions(-) diff --git a/Benchmark/NearestNeighbor/main.cpp b/Benchmark/NearestNeighbor/main.cpp index 319052fe..706a93cd 100644 --- a/Benchmark/NearestNeighbor/main.cpp +++ b/Benchmark/NearestNeighbor/main.cpp @@ -19,7 +19,9 @@ // halves that SIMD width -- the main reason it does not reach 2x in this table. It is also compiled // without -fopenmp, so its build is single-threaded like everything else here; kd3's docs quote a // much faster build, but that number is OpenMP-parallel (its author notes it is still faster serially). -// So these kd3 numbers are a fair single-threaded comparison, understated only by the double precision. +// So the `kd3 (double)` numbers are a fair single-threaded comparison, understated only by precision; +// `kd3 (float)` is additionally reported as kd3's native SoA/SIMD best case -- a different-precision +// reference (like TriangleMeshDistance's double column in the MeshSDF benchmark), not apples-to-apples. #include #include @@ -46,6 +48,9 @@ using Vec3 = EBGeometry::Vec3T; // kd3 tree, in double precision (its distance_t defaults to float; pin it to T for a fair comparison). using Kd3 = kd3::KdTree>; +// kd3 tree in its native float precision -- its SoA/SIMD best case (not a same-precision comparison). +using Kd3f = kd3::KdTree>; + constexpr std::size_t numPoints = 500000; constexpr std::size_t sampleSize = 500; constexpr std::uint64_t pointSeed = 123456789ULL; @@ -295,7 +300,7 @@ runCase(const std::string& a_label, const std::vector& a_positions) const double buildMs = 1.0e3 * timer.seconds(); if (!treeExpected) { - row("kd3", buildMs, std::numeric_limits::infinity(), sampleSize); // build failed -> flag it + row("kd3 (double)", buildMs, std::numeric_limits::infinity(), sampleSize); // build failed } else { const Kd3& tree = *treeExpected; @@ -327,7 +332,64 @@ runCase(const std::string& a_label, const std::vector& a_positions) for (std::size_t s = 0; s < sampleSize; s++) { bad += !ok(nnOther(s * stride), s); } - row("kd3", buildMs, queryUs, bad); + row("kd3 (double)", buildMs, queryUs, bad); + } + } + + // ── kd3 in native float precision (its SoA/SIMD best case; NOT a same-precision comparison) ── + { + std::vector fat(n); + timer.start(); + for (std::size_t i = 0; i < n; i++) { + fat[i] = Kd3f::FatPoint{{static_cast(a_positions[i][0]), + static_cast(a_positions[i][1]), + static_cast(a_positions[i][2])}, + static_cast(i)}; + } + + auto treeExpected = Kd3f::build(fat); + timer.stop(); + const double buildMs = 1.0e3 * timer.seconds(); + + if (!treeExpected) { + row("kd3 (float)", buildMs, std::numeric_limits::infinity(), sampleSize); // build failed + } + else { + const Kd3f& tree = *treeExpected; + + auto nnOther = [&](std::size_t i) -> float { + const Kd3f::point_t q = {static_cast(a_positions[i][0]), + static_cast(a_positions[i][1]), + static_cast(a_positions[i][2])}; + std::array buf{}; + const auto res = tree.query_knn(q, buf); + if (res) { + for (const auto& kr : *res) { + if (kr.payload_id != static_cast(i)) { + return kr.dist_sq; + } + } + } + return std::numeric_limits::max(); + }; + + volatile float sink = 0.0f; + timer.start(); + for (const std::uint32_t p : order) { + sink += nnOther(p); + } + timer.stop(); + (void)sink; + const double queryUs = 1.0e6 * timer.seconds() / double(n); + + // Float-precision cross-check: looser tolerance than the double methods (float has ~7 digits), + // since the returned squared distance is computed in float against the double brute-force truth. + std::size_t bad = 0; + for (std::size_t s = 0; s < sampleSize; s++) { + const double got = static_cast(nnOther(s * stride)); + bad += (std::abs(got - truth[s]) > 1.0e-4 * std::max(truth[s], 1.0)) ? 1 : 0; + } + row("kd3 (float)", buildMs, queryUs, bad); } } diff --git a/Benchmark/README.md b/Benchmark/README.md index 923197b9..5c096a83 100644 --- a/Benchmark/README.md +++ b/Benchmark/README.md @@ -54,7 +54,8 @@ PointCloudBVH ~118 0.62 ~105 0.48 PointCloudHashGrid ~12 1.50 ~11 1.60 picoflann ~104 0.75 ~104 0.48 nanoflann ~225 0.41 ~230 0.36 -kd3 ~83 0.57 ~86 0.50 +kd3 (double) ~83 0.57 ~82 0.45 +kd3 (float) ~75 0.51 ~74 0.42 ``` - The KD-tree queries are iterated in **Hilbert order** so their node cache is as warm as @@ -79,6 +80,12 @@ kd3 ~83 0.57 ~86 0.50 OpenMP-parallel (its author notes it is still faster serially). In this fair single-threaded mode kd3 posts the fastest build and a competitive query. (kd3 needs C++23, so the benchmark is built with `-std=c++23`.) +- **`kd3 (float)` is additionally shown as kd3's native best case** — its SoA/SIMD fast path in its + intended precision (a different-precision reference, not apples-to-apples with the double field, + exactly like TriangleMeshDistance's double column in the MeshSDF benchmark). It runs ~7–11% faster + than `kd3 (double)` here; the gain is modest on this particular machine because its AVX-512 clocks + down, damping the float SIMD-width advantage — elsewhere the float/double gap (and kd3's lead) is + larger. This is precisely why the tables are labelled machine-dependent snapshots. `MeshSDF/` — closest-point on a triangle mesh --------------------------------------------- diff --git a/Docs/Sphinx/source/Benchmark.rst b/Docs/Sphinx/source/Benchmark.rst index ba38bdb2..7be76f79 100644 --- a/Docs/Sphinx/source/Benchmark.rst +++ b/Docs/Sphinx/source/Benchmark.rst @@ -30,8 +30,8 @@ the full detail and representative numbers. * :file:`Benchmark/NearestNeighbor` -- all-nearest-neighbor over a point cloud: ``PointCloudBVH`` and ``PointCloudHashGrid`` vs `nanoflann `_ vs `picoflann `_ vs - `kd3 `_ (a SoA/SIMD KD-tree; needs C++23, run here in double -- - halving its float SIMD width -- for a same-precision comparison). + `kd3 `_ (a SoA/SIMD KD-tree; needs C++23, shown both in double + for a same-precision comparison and in its native float as a best-case reference). * :file:`Benchmark/MeshSDF` -- closest-point on a triangle mesh: ``TriMeshSDF`` vs `fcpw `_ (built with its Enoki CPU vectorization) vs `TriangleMeshDistance `_. From db8de238d845ddb9790e281c152554b72712871c Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 23 Jul 2026 15:27:40 +0000 Subject: [PATCH 17/17] Make disabled EBGEOMETRY_EXPECT truly vanish (unevaluated sizeof) The disabled-assertions branch expanded to ((void)(cond)), which still *evaluates* cond -- relying on the optimiser to dead-code-eliminate it. That elimination is not guaranteed: any assertion condition the compiler can't prove side-effect-free (e.g. the std::isfinite() precondition checks in the SoA distance kernels, evaluated on the order of 20x per nearest-neighbor query) stayed in the generated code. Switch the disabled expansion to (static_cast(sizeof((cond)))). sizeof is an unevaluated context, so cond is parsed -- keeping it syntax-checked and keeping assertion-only variables/parameters from tripping unused-entity warnings -- but is never executed, at any optimisation level. Disabled assertions now have exactly zero runtime cost and cannot have side effects. Verified: InstantiateAll (all classes, both precisions) compiles clean under -Wall -Wextra with assertions off; debug (assertions on) and release-test (assertions off) test suites both pass. An interleaved A/B of the NearestNeighbor benchmark shows PointCloudBVH's uniform-cube query ~4-5% faster with the expressions truly gone. Update ConfigurationOptions.rst to match. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01JhdRJydutVJqQjg2Fuee9Z --- Docs/Sphinx/source/ConfigurationOptions.rst | 13 ++++++++----- Source/EBGeometry_Macros.hpp | 15 ++++++++++----- 2 files changed, 18 insertions(+), 10 deletions(-) diff --git a/Docs/Sphinx/source/ConfigurationOptions.rst b/Docs/Sphinx/source/ConfigurationOptions.rst index 4f69e026..69715345 100644 --- a/Docs/Sphinx/source/ConfigurationOptions.rst +++ b/Docs/Sphinx/source/ConfigurationOptions.rst @@ -76,11 +76,14 @@ When ``EBGEOMETRY_ENABLE_ASSERTIONS`` is **not** defined (the default): .. code-block:: cpp - #define EBGEOMETRY_EXPECT(cond) ((void)(cond)) - -The condition is evaluated (preventing unused-variable warnings) but the branch is -absent from the generated code — a modern optimising compiler eliminates it entirely -at ``-O2`` or higher. + #define EBGEOMETRY_EXPECT(cond) (static_cast(sizeof((cond)))) + +The condition is **not** evaluated: ``sizeof`` is an unevaluated context, so the +expression is parsed (which keeps it syntax-checked, and keeps variables or +parameters that appear only inside assertions from tripping unused-entity warnings) +but is never executed. Disabled assertions therefore have exactly zero runtime cost +and cannot produce side effects, at any optimisation level — not merely once the +optimiser eliminates a discarded branch. When ``EBGEOMETRY_ENABLE_ASSERTIONS`` **is** defined: diff --git a/Source/EBGeometry_Macros.hpp b/Source/EBGeometry_Macros.hpp index c532c92e..064baed9 100644 --- a/Source/EBGeometry_Macros.hpp +++ b/Source/EBGeometry_Macros.hpp @@ -24,10 +24,13 @@ * if it is false, prints a diagnostic message to @c stderr and calls * @c std::abort(). * - * When @c EBGEOMETRY_ENABLE_ASSERTIONS is not defined the macro still - * evaluates the condition (to suppress "unused variable" warnings from - * variables that appear only inside assertions) but the result is - * discarded at zero runtime cost. + * When @c EBGEOMETRY_ENABLE_ASSERTIONS is not defined the macro does + * @b not evaluate the condition at all: it expands to a discarded + * @c sizeof, which is an unevaluated context. The expression is therefore + * never executed (guaranteed zero runtime cost and no side effects, even + * at @c -O0), yet is still parsed -- so it stays syntax-checked and any + * variables/parameters that appear only inside assertions are still + * considered "used" and do not trigger unused-entity warnings. * * @par Enabling assertions * @code{.cmake} @@ -59,7 +62,9 @@ } \ } while (0) #else -#define EBGEOMETRY_EXPECT(cond) ((void)(cond)) +// Unevaluated sizeof: cond is parsed (syntax-checked, names count as "used") but never evaluated, +// so disabled assertions have exactly zero runtime cost and cannot have observable side effects. +#define EBGEOMETRY_EXPECT(cond) (static_cast(sizeof((cond)))) #endif #endif // EBGEOMETRY_MACROS_HPP