diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index e4291c2e..7891cda9 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 Benchmark/kd3 + - 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 1915074c..486b857a 100644 --- a/.gitmodules +++ b/.gitmodules @@ -2,3 +2,22 @@ path = common-3d-test-models url = https://github.com/alecjacobson/common-3d-test-models shallow = true +[submodule "Benchmark/nanoflann"] + path = Benchmark/nanoflann + url = https://github.com/jlblancoc/nanoflann.git + shallow = true +[submodule "Benchmark/picoflann"] + path = Benchmark/picoflann + url = https://github.com/rmsalinas/picoflann.git + shallow = true +[submodule "Benchmark/fcpw"] + path = Benchmark/fcpw + url = https://github.com/rohan-sawhney/fcpw.git + shallow = true +[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/MeshSDF/GNUmakefile b/Benchmark/MeshSDF/GNUmakefile new file mode 100644 index 00000000..b6150a1c --- /dev/null +++ b/Benchmark/MeshSDF/GNUmakefile @@ -0,0 +1,32 @@ +# 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 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. + +EBGEOMETRY_HOME ?= ../.. +CXX ?= g++ +CXXFLAGS ?= -std=c++17 -O3 -march=native +FCPW_SIMD_WIDTH ?= 8 + +INCLUDES := -I$(EBGEOMETRY_HOME) \ + -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) + +TARGET := MeshSDF.ex + +$(TARGET): main.cpp + $(CXX) $(CXXFLAGS) $(DEFINES) $(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..479ea5ef --- /dev/null +++ b/Benchmark/MeshSDF/main.cpp @@ -0,0 +1,182 @@ +// SPDX-FileCopyrightText: 2026 Robert Marskar +// +// SPDX-License-Identifier: GPL-3.0-or-later + +// 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 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 +#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] : "../../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, true /* vectorize (Enoki MBVH) */); + 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); + + // ── 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))) { + badFcpw++; + } + if (std::abs(ebDist[i] - tmdDist[i]) > T(1.0e-3) * std::max(tmdDist[i], T(1))) { + badTmd++; + } + } + + 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/NearestNeighbor/GNUmakefile b/Benchmark/NearestNeighbor/GNUmakefile new file mode 100644 index 00000000..54ecf841 --- /dev/null +++ b/Benchmark/NearestNeighbor/GNUmakefile @@ -0,0 +1,27 @@ +# 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++23 -O3 -march=native + +INCLUDES := -I$(EBGEOMETRY_HOME) \ + -I$(EBGEOMETRY_HOME)/Benchmark/picoflann \ + -I$(EBGEOMETRY_HOME)/Benchmark/nanoflann/include \ + -I$(EBGEOMETRY_HOME)/Benchmark/kd3/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..706a93cd --- /dev/null +++ b/Benchmark/NearestNeighbor/main.cpp @@ -0,0 +1,412 @@ +// SPDX-FileCopyrightText: 2026 Robert Marskar +// +// SPDX-License-Identifier: GPL-3.0-or-later + +// 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. 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 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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>; + +// 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; + +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>; + +// 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) +{ + 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; +} + +// 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) +{ + 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(a_positions); + timer.stop(); + const double sortUsPerPt = 1.0e6 * timer.seconds() / double(n); + + 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, a_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(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(n); + + 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); + } + + // ── 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(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(a_positions, a_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(n); + + std::size_t bad = 0; + for (std::size_t s = 0; s < sampleSize; s++) { + 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) { + got = pr.second; + break; + } + } + bad += !ok(got, s); + } + row("picoflann", buildMs, queryUs, bad); + } + + // ── nanoflann (per-point findNeighbors, k=2) ── + { + NanoCloud cloud{a_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] = {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]; + }; + + 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("nanoflann", buildMs, queryUs, bad); + } + + // ── 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 (double)", buildMs, std::numeric_limits::infinity(), sampleSize); // build failed + } + 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 (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); + } + } + + 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"; +} + +} // namespace + +int +main() +{ + 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)); + + return 0; +} diff --git a/Benchmark/README.md b/Benchmark/README.md new file mode 100644 index 00000000..5c096a83 --- /dev/null +++ b/Benchmark/README.md @@ -0,0 +1,114 @@ +Benchmark +========= + +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 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) +* [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/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) +``` + +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 — correctness is the part +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 cloud (double precision), find its nearest *other* point. Compares +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: + +* **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): + +``` + uniform cube sphere surface +Method Build(ms) Query(us/pt) Build(ms) Query(us/pt) +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 (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 + 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 ~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'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`.) +- **`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 +--------------------------------------------- + +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 ~63 2.6 +fcpw (Enoki) ~46 3.5 +TriangleMeshDistance ~115 10.5 +``` + +- **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/Benchmark/TriangleMeshDistance b/Benchmark/TriangleMeshDistance new file mode 160000 index 00000000..5530eaa8 --- /dev/null +++ b/Benchmark/TriangleMeshDistance @@ -0,0 +1 @@ +Subproject commit 5530eaa85537cc62b884c3f8a3e64e000129eb93 diff --git a/Benchmark/fcpw b/Benchmark/fcpw new file mode 160000 index 00000000..61814ff0 --- /dev/null +++ b/Benchmark/fcpw @@ -0,0 +1 @@ +Subproject commit 61814ff0c7b69d61dac3b9725cda1541b1b3ec4f 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/Benchmark/nanoflann b/Benchmark/nanoflann new file mode 160000 index 00000000..ff0eb50a --- /dev/null +++ b/Benchmark/nanoflann @@ -0,0 +1 @@ +Subproject commit ff0eb50a4b972c642c54fe453095b724d9ae7c60 diff --git a/Benchmark/picoflann b/Benchmark/picoflann new file mode 160000 index 00000000..d5fd165a --- /dev/null +++ b/Benchmark/picoflann @@ -0,0 +1 @@ +Subproject commit d5fd165aa8f0f91ceb6277252db3e8f4455cb470 diff --git a/Docs/Sphinx/source/Benchmark.rst b/Docs/Sphinx/source/Benchmark.rst new file mode 100644 index 00000000..7be76f79 --- /dev/null +++ b/Docs/Sphinx/source/Benchmark.rst @@ -0,0 +1,37 @@ +.. _Chap:Benchmark: + +Benchmarks +========== + +.. important:: + + 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 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`` and + ``PointCloudHashGrid`` vs `nanoflann `_ vs + `picoflann `_ vs + `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 `_. 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/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/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 99f62e58..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 at the repository root (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 @@ -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 = [ 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