forked from root-project/root
-
Notifications
You must be signed in to change notification settings - Fork 10
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
[geom] New BVH/VoxelGrid-based implementation of TGeoParallelWorld ro…
…utines This commit provides a rewrite of key functions of TGeoParallelWorld, achieving: (a) faster initialization (b) faster execution of the Safety function (from ~O(N) to ~O(1)) (c) similar execution of FindNode/FindBoundary functions (~log(N)) (d) less memory consumption (better memory scalability) The development for this commit was motivated from a use case in ALICE, in which the parallel world "scene" can be very large (~100K volumes). In this case, TGeoVoxelFinder takes very long to construct and consumes a very large amount of memory (GBs). In addition, the evaluation of the Safety function dominates the Geant simulation time. The improvements in this commit are mainly achieved through: * The use of a boundary volume hierarchy (BVH) as the base acceleration entity, replacing TGeoVoxelFinder. BVH are the standard in industry/computer-graphics, for what concerns ray-object intersection tasks. The BVH is constructed from axis-aligned bounding boxes and employed in the FindBoundary/FindNode implementations * The use of a 3D voxel grid (TGeoVoxelGrid) structure, able to store properties "local" or in the vicinity of a cartesian coordinate P. This structure allows to reduce the (typical) algorithmic complexity for "Safety" queries to ~O(1) (with a constant factor determined by the voxel size). Filling of the 3D voxel grid cache for Safety is done on-the-fly (using the BVH once). Ideas for these improvements come from prior work in related libraries such as VecGeom. ----- Implementation details: * The implementation is, for now (until fully tested), provided in a backward compatible manner: - By default, nothing changes - Users have to activate the BVH mode via TGeoParallelWorld::SetAccelerationMode - Users may hence compare the 2 modes or choose the best depending on complexity and needs * Functions for Safety, FindNode, FindBoundary dispatch to some internal implementation. This causes an extra lookup/jump, which can be removed once BVH is fully validated * For the BVH, a well known open source implementation is included in header-only form. The headers are copied from https://github.com/madmann91/bvh commit 66e445b92f68801a6dd8ef943fe3038976ecb4ff. Some minor patching has been done in order to achieve compilation for C++17 (see README). * A new class, TGeoVoxelGrid is provided for the cartesian VoxelGrid container. ---- Performance examples: In a test with the ALICE simulation framework including the ITS + TPC detectors with 48240 volumes on the parallel world, we see * initialization time goes from TGeoVoxelFinder: 10s ---> BVH: 40ms * Geant simulation time: 10s --> 2s * memory usage: 3GB --> 1GB Hence, this PR will make a big difference for the ALICE simulation program. It was verified, that identical results (number of hits, steps, etc) are obtained when going from TGeoVoxelFinder --> BVH+GRID. ---- Outlook: Similar techniques could be applied to ordinary TGeoNavigator routines.
- Loading branch information
Showing
24 changed files
with
3,748 additions
and
10 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,176 @@ | ||
| /// \author Sandro Wenzel <sandro.wenzel@cern.ch> | ||
| /// \date 2024-02-22 | ||
|
|
||
| /************************************************************************* | ||
| * Copyright (C) 1995-2024, Rene Brun and Fons Rademakers. * | ||
| * All rights reserved. * | ||
| * * | ||
| * For the licensing terms see $ROOTSYS/LICENSE. * | ||
| * For the list of contributors see $ROOTSYS/README/CREDITS. * | ||
| *************************************************************************/ | ||
|
|
||
| #ifndef ROOT_TGeoVoxelGrid | ||
| #define ROOT_TGeoVoxelGrid | ||
|
|
||
| // a simple structure to encode voxel indices, to address | ||
| // individual voxels in the 3D grid. | ||
| struct TGeoVoxelGridIndex { | ||
| int ix{-1}; | ||
| int iy{-1}; | ||
| int iz{-1}; | ||
| size_t idx{std::numeric_limits<size_t>::max()}; | ||
| bool isValid() const { return idx != std::numeric_limits<size_t>::max(); } | ||
| }; | ||
|
|
||
| /// A finite 3D grid structure, mapping/binning arbitrary 3D cartesian points | ||
| /// onto discrete "voxels". Each such voxel can store an object of type T. | ||
| /// The precision of the voxel binning is done with S (float or double). | ||
| template <typename T, typename S = float> | ||
| class TGeoVoxelGrid { | ||
| public: | ||
| TGeoVoxelGrid(S xmin, S ymin, S zmin, S xmax, S ymax, S zmax, S Lx_, S Ly_, S Lz_) | ||
| : fMinBound{xmin, ymin, zmin}, fMaxBound{xmax, ymax, zmax}, fLx(Lx_), fLy(Ly_), fLz(Lz_) | ||
| { | ||
|
|
||
| // Calculate the number of voxels in each dimension | ||
| fNx = static_cast<int>((fMaxBound[0] - fMinBound[0]) / fLx); | ||
| fNy = static_cast<int>((fMaxBound[1] - fMinBound[1]) / fLy); | ||
| fNz = static_cast<int>((fMaxBound[2] - fMinBound[2]) / fLz); | ||
|
|
||
| finvLx = 1. / fLx; | ||
| finvLy = 1. / fLy; | ||
| finvLz = 1. / fLz; | ||
|
|
||
| fHalfDiag = std::sqrt(fLx / 2. * fLx / 2. + fLy / 2. * fLy / 2. + fLz / 2. * fLz / 2.); | ||
|
|
||
| // Resize the grid to hold the voxels | ||
| fGrid.resize(fNx * fNy * fNz); | ||
| } | ||
|
|
||
| T &at(int i, int j, int k) { return fGrid[index(i, j, k)]; } | ||
|
|
||
| // check if point is covered by voxel structure | ||
| bool inside(std::array<S, 3> const &p) const | ||
| { | ||
| for (int i = 0; i < 3; ++i) { | ||
| if (p[i] < fMinBound[i] || p[i] > fMaxBound[i]) { | ||
| return false; | ||
| } | ||
| } | ||
| return true; | ||
| } | ||
|
|
||
| // Access a voxel given a 3D point P | ||
| T &at(std::array<S, 3> const &P) | ||
| { | ||
| int i, j, k; | ||
| pointToVoxelIndex(P, i, j, k); // Convert point to voxel index | ||
| return fGrid[index(i, j, k)]; // Return reference to voxel's data | ||
| } | ||
|
|
||
| T *at(TGeoVoxelGridIndex const &vi) | ||
| { | ||
| if (!vi.isValid()) { | ||
| return nullptr; | ||
| } | ||
| return &fGrid[vi.idx]; | ||
| } | ||
|
|
||
| // Set the data of a voxel at point P | ||
| void set(std::array<S, 3> const &p, const T &value) | ||
| { | ||
| int i, j, k; | ||
| pointToVoxelIndex(p, i, j, k); // Convert point to voxel index | ||
| fGrid[index(i, j, k)] = value; // Set the value at the voxel | ||
| } | ||
|
|
||
| // Set the data of a voxel at point P | ||
| void set(int i, int j, int k, const T &value) | ||
| { | ||
| fGrid[index(i, j, k)] = value; // Set the value at the voxel | ||
| } | ||
|
|
||
| void set(TGeoVoxelGridIndex const &vi, const T &value) { fGrid[vi.idx] = value; } | ||
|
|
||
| // Get voxel dimensions | ||
| int getVoxelCountX() const { return fNx; } | ||
| int getVoxelCountY() const { return fNy; } | ||
| int getVoxelCountZ() const { return fNz; } | ||
|
|
||
| // returns the cartesian mid-point coordinates of a voxel given by a VoxelIndex | ||
| std::array<S, 3> getVoxelMidpoint(TGeoVoxelGridIndex const &vi) const | ||
| { | ||
| const S midX = fMinBound[0] + (vi.ix + 0.5) * fLx; | ||
| const S midY = fMinBound[1] + (vi.iy + 0.5) * fLy; | ||
| const S midZ = fMinBound[2] + (vi.iz + 0.5) * fLz; | ||
|
|
||
| return {midX, midY, midZ}; | ||
| } | ||
|
|
||
| S getDiagonalLength() const { return fHalfDiag; } | ||
|
|
||
| // Convert a point p(x, y, z) to voxel indices (i, j, k) | ||
| // if point is outside set indices i,j,k to -1 | ||
| void pointToVoxelIndex(std::array<S, 3> const &p, int &i, int &j, int &k) const | ||
| { | ||
| if (!inside(p)) { | ||
| i = -1; | ||
| j = -1; | ||
| k = -1; | ||
| } | ||
|
|
||
| i = static_cast<int>((p[0] - fMinBound[0]) * finvLx); | ||
| j = static_cast<int>((p[1] - fMinBound[1]) * finvLy); | ||
| k = static_cast<int>((p[2] - fMinBound[2]) * finvLz); | ||
|
|
||
| // Clamp the indices to valid ranges | ||
| i = std::min(i, fNx - 1); | ||
| j = std::min(j, fNy - 1); | ||
| k = std::min(k, fNz - 1); | ||
| } | ||
|
|
||
| // Convert a point p(x, y, z) to voxel index object | ||
| // if outside, an invalid index object will be returned | ||
| TGeoVoxelGridIndex pointToVoxelIndex(std::array<S, 3> const &p) const | ||
| { | ||
| if (!inside(p)) { | ||
| return TGeoVoxelGridIndex(); // invalid voxel index | ||
| } | ||
|
|
||
| int i = static_cast<int>((p[0] - fMinBound[0]) * finvLx); | ||
| int j = static_cast<int>((p[1] - fMinBound[1]) * finvLy); | ||
| int k = static_cast<int>((p[2] - fMinBound[2]) * finvLz); | ||
|
|
||
| // Clamp the indices to valid ranges | ||
| i = std::min(i, fNz - 1); | ||
| j = std::min(j, fNy - 1); | ||
| k = std::min(k, fNz - 1); | ||
|
|
||
| return TGeoVoxelGridIndex{i, j, k, index(i, j, k)}; | ||
| } | ||
|
|
||
| TGeoVoxelGridIndex pointToVoxelIndex(S x, S y, S z) const { return pointToVoxelIndex(std::array<S, 3>{x, y, z}); } | ||
|
|
||
| // Convert voxel indices (i, j, k) to a linear index in the grid array | ||
| size_t index(int i, int j, int k) const { return i + fNx * (j + fNy * k); } | ||
|
|
||
| void indexToIndices(size_t idx, int &i, int &j, int &k) const | ||
| { | ||
| k = idx % fNz; | ||
| j = (idx / fNz) % fNy; | ||
| i = idx / (fNy * fNz); | ||
| } | ||
|
|
||
| // member data | ||
|
|
||
| std::array<S, 3> fMinBound; | ||
| std::array<S, 3> fMaxBound; // 3D bounds for grid structure | ||
| S fLx, fLy, fLz; // Voxel dimensions | ||
| S finvLx, finvLy, finvLz; // inverse voxel dimensions | ||
| S fHalfDiag; // cached value for voxel half diagonal length | ||
|
|
||
| int fNx, fNy, fNz; // Number of voxels in each dimension | ||
| std::vector<T> fGrid; // The actual voxel grid data container | ||
| }; | ||
|
|
||
| #endif |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,23 @@ | ||
| add_library(bvh INTERFACE) | ||
| find_package(Threads) | ||
| if (Threads_FOUND) | ||
| # Link with the threading library of the system, which may | ||
| # be required by standard header <thread> on some systems | ||
| target_link_libraries(bvh INTERFACE Threads::Threads) | ||
| endif() | ||
|
|
||
| target_include_directories(bvh INTERFACE | ||
| $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/src> | ||
| $<INSTALL_INTERFACE:include>) | ||
|
|
||
| set_target_properties(bvh PROPERTIES CXX_STANDARD 20) | ||
|
|
||
| install( | ||
| DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/ | ||
| DESTINATION include/bvh/v2 | ||
| FILES_MATCHING PATTERN "*.h" | ||
| PATTERN "c_api" EXCLUDE) | ||
|
|
||
| if (BVH_BUILD_C_API) | ||
| add_subdirectory(c_api) | ||
| endif() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,13 @@ | ||
| These headers, providing routines to construct and navigate | ||
| bounding volume hierarchies, have been copied from https://github.com/madmann91/bvh | ||
| commit 66e445b92f68801a6dd8ef94. | ||
|
|
||
|
|
||
| Minor changes have been subsequently been applied to achieve compilation with C++17: | ||
|
|
||
| - inclusion of alternative span, when std::span is not found | ||
| - replacement of C++20 defaulted comparison operators with actual implementation | ||
| - old-style struct construction for objects of type "Reinsertion" | ||
| - use of std::inner_product instead of std::transform_reduce (gcc 8.5 had problems) | ||
|
|
||
| This is needed since ROOT should compile with C++17. |
Oops, something went wrong.