Add BVH-based navigation functions to TGeoTessellated #21045

sawenzel · 2026-01-27T13:12:42Z

This commit adds implementations of the navigation interfaces required by TGeoNavigator to TGeoTessellated.

As a result, TGeoTessellated can now be used in applications that rely on this navigator, such as VMC-based GEANT simulations.

The implementation is based on fast Bounding Volume Hierarchy (BVH) queries introduced in a recent TGeoParallelWorld rewrite.

The BVH is only exposed to the internal implementation. It is not streamed to disk for now and therefore needs to be reinitialized when reading from TBuffer. For this reason, a custom streamer for TGeoTessellated is proposed.

The implementation shows very good performance, being multiple times faster than comparable implementations in Geant4 or VecGeom. The latter were also used to validate the algorithms. Navigation functionality has additionally been tested at a high level using VecGeom’s XRayNavigationBenchmarker.

Note that this is an initial implementation and should be considered beta quality. The following limitations and remarks apply:

DistFromInside and DistFromOutside do not implement all iact= cases. Only iact == 3 is supported, which is the case actually used by TGeoNavigator.
The algorithms have not yet been validated in large-scale production environments.
Further optimization seem possible and will be done in future commits.

This work is motivated by the need to GEANT-simulate detector prototypes from CAD models in ALICE. ALICE is using TGeo.

Finally, the commit also fixes some minor bugs in the TGeoChecker with respect to bounding box origins, previously not taken into account.

This commit adds implementations of the navigation interfaces required by TGeoNavigator to TGeoTessellated. As a result, TGeoTessellated can now be used in applications that rely on this navigator, such as VMC-based GEANT simulations. The implementation is based on fast Bounding Volume Hierarchy (BVH) queries introduced in a recent TGeoParallelWorld rewrite. The BVH is only exposed to the internal implementation. It is not streamed to disk for now and therefore needs to be reinitialized when reading from TBuffer. For this reason, a custom streamer for TGeoTessellated is proposed. The implementation shows very good performance, being multiple times faster than comparable implementations in Geant4 or VecGeom. The latter were also used to validate the algorithms. Navigation functionality has additionally been tested at a high level using VecGeom’s XRayNavigationBenchmarker. Note that this is an initial implementation and should be considered beta quality. The following limitations and remarks apply: - DistFromInside and DistFromOutside do not implement all `iact=` cases. Only `iact == 3` is supported, which is the case actually used by TGeoNavigator. - The algorithms have not yet been validated in large-scale production environments. - Further optimization seem possible and will be done in future commits. This work is motivated by the need to GEANT-simulate detector prototypes from CAD models in ALICE. ALICE is using TGeo. Finally, the commit also fixes some minor bugs in the TGeoChecker with respect to bounding box origins, preveously not taken into account.

sawenzel · 2026-01-27T13:16:13Z

Sharing code for initial discussion. Adding a table to get idea about performance for a 20K triangular faceted sphere (testobject):

The simple BVH performs well in comparison to the Voxelization implementation of G4 and the SIMD-based structure of VecGeom.

sawenzel · 2026-01-27T13:20:29Z

Further, validation tests have been performed in the level of the global navigation system (G4Navigator, TGeoNavigator, VecGeomNavigator) within the XRayBenchmark of VecGeom. All implementations yield the same XRay image which gives confidence that the present implementation is ok.

Overall, the present PR should present a significant step forward for the TGeo navigation system to handle tessellated shapes from CAD models.

geom/geom/inc/LinkDef1.h

geom/geom/inc/TGeoTessellated.h

dpiparo · 2026-01-28T08:24:05Z

geom/geom/src/TGeoTessellated.cxx


-   if (fNfacets > 0 && GetNfacets() == fNfacets)
-      CloseShape();
+   // if (fNfacets > 0 && GetNfacets() == fNfacets)


I would remove the dead code if possible

dpiparo · 2026-01-28T08:25:06Z

geom/geom/src/TGeoTessellated.cxx

+   return (t > rayEPS) ? t : INF;
+}
+
+template <typename T = float>


Would ROOT::Math::XYZVector be usable here instead?

It could, if we replace the typedef Tessellated::Vertex_t point to ROOT::Math::XYZVector, with probably some tweaks, because the interfaces are not identical. Probably this is desirable for future maintainability, but IMO can/should be done in a separate PR. And XYZVector could benefit from some simple additions, such as magnitude (now only squared available), or indexed access.

dpiparo · 2026-01-28T08:26:56Z

geom/geom/src/TGeoTessellated.cxx

+
+   // we need bounding boxes enclosing the primitives and centers of primitives
+   // (replaced here by centers of bounding boxes) to build the bvh
+   auto bboxes_ptr = new std::vector<BBox>();


I do not understand this allocation. Why? Who deletes this memory?

It is assigned to fBVH, but not deleted in the dtor

dpiparo · 2026-01-28T08:27:31Z

geom/geom/src/TGeoTessellated.cxx

+   // we need bounding boxes enclosing the primitives and centers of primitives
+   // (replaced here by centers of bounding boxes) to build the bvh
+   auto bboxes_ptr = new std::vector<BBox>();
+   // fBoundingBoxes = (void *)bboxes_ptr;


If we don't need this line, perhaps we can remove it instead of commenting?

dpiparo · 2026-01-28T08:28:59Z

geom/geom/src/TGeoTessellated.cxx

+   }
+
+   // doesn't need to be normalized and probes all normals
+   double test_dir[3] = {1.0, 1.41421356237, 1.73205080757};


maybe a comment about these numbers is in order?

Right, we should strive to document better the implementation so that it will be easier to maintain. To start with, we can do it in a doxygen @details banner for the methods containing sensibly-complex algorithms

dpiparo · 2026-01-28T08:30:28Z

geom/geom/src/TGeoTessellated.cxx

+      vol +=
+         a[0] * (b[1] * c[2] - b[2] * c[1]) + b[0] * (c[1] * a[2] - c[2] * a[1]) + c[0] * (a[1] * b[2] - a[2] * b[1]);
+   }
+   return vol / 6.0;


oppportunity to get rid of a division?

dpiparo · 2026-01-28T08:34:49Z

geom/geom/src/TGeoTessellated.cxx

+// A priority queue for BVHPrioElement with an additional clear method
+// for quick reset
+template <typename Comparator>
+class BVHPrioQueue : public std::priority_queue<BVHPrioElement, std::vector<BVHPrioElement>, Comparator> {


Can we find an alternative to subclassing stl containers?

dpiparo · 2026-01-28T08:39:45Z

Thanks a lot for this PR. I left a few comments. In general, I would like to ask whether you think it's possible to add some test coverage for this new nice functionality. The coverage in the geometry package in general can be improved, and here we have an opportunity!

agheata

Thanks for this great development and congrats for the results @sawenzel ! I have made some comments inline. Here are a few more general comments - also as a personal reminder for the needed follow-ups - most of them do not need to be addressed in this PR, but are important.

Persistency might become important for the optimization structure. I don't know if you have benchmarked the scaling of the bvh creation time with the number of triangles. IMO we should make bvh::v2 persistence-friendly and remove the custom streamer this. Not necessary in this PR, but very soon, because adopting a persistent layout tends to be forever.
We should strive to reuse the vector types. Legacy TGeo is double[3] (terrible), then I introduced TGeoVector3, aliased as Vertex_t to have a decent interface, ignoring that ROOT comes with XYZVector, then you added bvh::v2 coming with its own vectors, and there are also some vector helpers added in this PR. We should make Vertex_t and bvh::v2::vec3 point to XYZVector, and possibly enhance the functionality of the latter with indexed access and a few other needed helper functions
The algorithms here (as well as in bvh::v2) are really missing documentation, which is quite a problem in TGeo anyway, so we should try to make the complex algorithms understandable. Extended doxygen per function is OK, I guess.
We should add unit & regression testing to make sure we have correctness and do not break it in the future. I am thinking of a ShapeTester utility in the VecGeom style, but there are other ways, since we can already replace ROOT solids with VecGeom ones and compare some output - even if this requires a version of ROOT built with VecGeom support. We could add something simple in this PR, such as a TestTessellated unit test comparing expected values for the different navigation functions, and also moving the _loop comparisons there. See tests in /geom/test
Longer term strategy, ROOT should strive to use VecGeom for navigation, so we will need this kind of implementation to be maintained there rather than in potentially three places.

Anyway, thanks for the great contribution!

agheata · 2026-02-02T08:31:25Z

geom/geom/inc/TGeoTessellated.h

+   Double_t DistFromInside_Loop(const Double_t *point, const Double_t *dir) const;
+   Double_t DistFromOutside_Loop(const Double_t *point, const Double_t *dir) const;
+   bool Contains_Loop(const Double_t *point) const;


Since these are for testing purposes only, leaving them in the public interface may not be the best choice. The alternative could be to move these to a test utility, where we can add also some regression tests in future. What about Safety_Loop ?

agheata · 2026-02-02T08:41:05Z

geom/geom/src/TGeoTessellated.cxx

 }
+
+// implementation of some geometry helper functions in anonymous namespace
+namespace {


Tessellated::Vertex_t is eventually TGeoVector3`, which already has these functions, so no need to re-implement them here

agheata · 2026-02-02T08:53:33Z

geom/geom/src/TGeoTessellated.cxx

+   return (t > rayEPS) ? t : INF;
+}
+
+template <typename T = float>


It could, if we replace the typedef Tessellated::Vertex_t point to ROOT::Math::XYZVector, with probably some tweaks, because the interfaces are not identical. Probably this is desirable for future maintainability, but IMO can/should be done in a separate PR. And XYZVector could benefit from some simple additions, such as magnitude (now only squared available), or indexed access.

agheata · 2026-02-02T08:55:17Z

geom/geom/src/TGeoTessellated.cxx

+}
+
+template <typename T = float>
+struct Vec3f {


Why this helper, it duplicates stuff, can we use Vertex_t?

agheata · 2026-02-02T09:01:19Z

geom/geom/src/TGeoTessellated.cxx

+   T invLen = T(1.0f) / std::sqrt(len2);
+   return {v.x * invLen, v.y * invLen, v.z * invLen};
+}
+


agheata · 2026-02-02T09:37:00Z

geom/geom/src/TGeoTessellated.cxx

+
+   // we need bounding boxes enclosing the primitives and centers of primitives
+   // (replaced here by centers of bounding boxes) to build the bvh
+   auto bboxes_ptr = new std::vector<BBox>();


It is assigned to fBVH, but not deleted in the dtor

agheata · 2026-02-02T09:41:27Z

geom/geom/src/TGeoTessellated.cxx

+   }
+
+   // doesn't need to be normalized and probes all normals
+   double test_dir[3] = {1.0, 1.41421356237, 1.73205080757};


Right, we should strive to document better the implementation so that it will be easier to maintain. To start with, we can do it in a doxygen @details banner for the methods containing sensibly-complex algorithms

agheata · 2026-02-02T09:46:25Z

geom/geom/src/TGeoTessellated.cxx

+   double local_step = Big();
+   // the ray used for bvh interaction
+   Ray ray(Vec3(point[0], point[1], point[2]),          // origin
+           Vec3(test_dir[0], test_dir[1], test_dir[2]), // direction


Is this normalized anywhere?

agheata · 2026-02-02T09:48:24Z

geom/geom/src/TGeoTessellated.cxx

+template <typename T>
+bool contains(bvh::v2::BBox<T, 3> const &box, bvh::v2::Vec<T, 3> const &p)
+{
+   auto min = box.min;
+   auto max = box.max;
+   return (p[0] >= min[0] && p[0] <= max[0]) && (p[1] >= min[1] && p[1] <= max[1]) &&
+          (p[2] >= min[2] && p[2] <= max[2]);
+}


It makes sense to put this in bvh::v2?

agheata · 2026-02-02T09:53:01Z

geom/geom/src/TGeoTessellated.cxx

+Double_t TGeoTessellated::DistFromInside_Loop(const Double_t *point, const Double_t *dir) const
+{
+   // Bias the starting point slightly along the direction
+   Vertex_t p(point[0], point[1], point[2]);
+   Vertex_t d(dir[0], dir[1], dir[2]);
+
+   double dist = Big();
+   for (size_t i = 0; i < fFacets.size(); ++i) {
+      const auto &facet = fFacets[i];
+      const auto &n = fOutwardNormals[i];
+
+      // Only exiting surfaces are relevant (from inside--> dot product must be positive)
+      if (n[0] * dir[0] + n[1] * dir[1] + n[2] * dir[2] <= 0.0) {
+         continue;
+      }
+
+      const auto &v0 = fVertices[facet[0]];
+      const auto &v1 = fVertices[facet[1]];
+      const auto &v2 = fVertices[facet[2]];
+
+      const double t = rayTriangle(p, d, v0, v1, v2, 0.);
+
+      if (t < dist) {
+         dist = t;
+      }
+   }
+   return dist;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// trivial (non-BVH) DistFromOutside function
+
+Double_t TGeoTessellated::DistFromOutside_Loop(const Double_t *point, const Double_t *dir) const
+{
+   // Bias the starting point slightly along the direction
+   Vertex_t p(point[0], point[1], point[2]);
+   Vertex_t d(dir[0], dir[1], dir[2]);
+
+   double dist = Big();
+   for (size_t i = 0; i < fFacets.size(); ++i) {
+      const auto &facet = fFacets[i];
+      const auto &n = fOutwardNormals[i];
+
+      // Only exiting surfaces are relevant (from outside, the dot product must be negative)
+      if (n[0] * dir[0] + n[1] * dir[1] + n[2] * dir[2] > 0.0) {
+         continue;
+      }
+
+      const auto &v0 = fVertices[facet[0]];
+      const auto &v1 = fVertices[facet[1]];
+      const auto &v2 = fVertices[facet[2]];
+
+      const double t = rayTriangle(p, d, v0, v1, v2, 0.);
+
+      if (t < dist) {
+         dist = t;
+      }
+   }
+   return dist;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// trivial (non-BVH) Contains
+
+bool TGeoTessellated::Contains_Loop(const Double_t *point) const
+{
+   // Fixed ray direction
+   const Vertex_t test_dir{1.0, 1.41421356237, 1.73205080757};
+
+   // Bias point slightly along ray to avoid t == 0 hits
+   Vertex_t p(point[0], point[1], point[2]);
+
+   int crossings = 0;
+   for (size_t i = 0; i < fFacets.size(); ++i) {
+      const auto &facet = fFacets[i];
+
+      const auto &v0 = fVertices[facet[0]];
+      const auto &v1 = fVertices[facet[1]];
+      const auto &v2 = fVertices[facet[2]];
+
+      const double t = rayTriangle(p, test_dir, v0, v1, v2, 0.);
+      if (t != std::numeric_limits<double>::infinity()) {
+         ++crossings;
+      }
+   }
+   return (crossings & 1);
+}


Candidates for a test utility, no need to keep these here

sawenzel · 2026-02-02T13:38:32Z

@dpiparo @agheata :Thanks for initial review. I'll take a look and see what can be addressed with reasonable effort.

geom/geom/src/TGeoTessellated.cxx

github-actions · 2026-02-02T17:29:00Z

Test Results

0 tests 0 ✅ 0s ⏱️
0 suites 0 💤
0 files 0 ❌

Results for commit 30199d6.

♻️ This comment has been updated with latest results.

GDML import was not building voxels because of early return if the check is disabled

silverweed · 2026-02-04T13:26:39Z

geom/geom/inc/TGeoTessellated.h

+
   std::multimap<long, int> fVerticesMap; //! Temporary map used to deduplicate vertices
+   bool fIsClosed = false;                //! to know if shape still needs closure/initialization
+   void *fBVH = nullptr;                  //! BVH acceleration structure for safety and navigation


If we put the real type, can we also use a unique_ptr while we're at it?

silverweed · 2026-02-04T13:28:34Z

geom/geom/src/TGeoTessellated.cxx

+
+   // loop over all the triangles/Facets;
+   int nd = fFacets.size();
+   for (int i = 0; i < nd; ++i) {


Since we know how many elements we're gonna push in:

Suggested change

for (int i = 0; i < nd; ++i) {

centers.reserve(nd);

bboxes.reserve(nd);

for (int i = 0; i < nd; ++i) {

silverweed · 2026-02-04T13:29:20Z

geom/geom/src/TGeoTessellated.cxx

+   *bvhptr = std::move(bvh); // copy structure
+   fBVH = (void *)(bvhptr);
+
+   return;


Suggested change

return;

sawenzel requested a review from agheata as a code owner January 27, 2026 13:12

dpiparo assigned agheata Jan 28, 2026

dpiparo reviewed Jan 28, 2026

View reviewed changes

geom/geom/inc/LinkDef1.h Show resolved Hide resolved

dpiparo reviewed Jan 28, 2026

View reviewed changes

geom/geom/inc/TGeoTessellated.h Outdated Show resolved Hide resolved

dpiparo reviewed Jan 28, 2026

View reviewed changes

BSali mentioned this pull request Feb 1, 2026

Implementation of the TGeoTessellated navigation functionality #21095

Draft

agheata requested changes Feb 2, 2026

View reviewed changes

agheata reviewed Feb 2, 2026

View reviewed changes

geom/geom/src/TGeoTessellated.cxx Outdated Show resolved Hide resolved

agheata and others added 7 commits February 3, 2026 14:54

[geom] Fix warnings in bvh/v2 headers.

dc87a3d

[geom] Fix warnings in TGeoTessellated

98d5515

[geom] Fix CloseShape handling of the check flag.

ebd09b6

GDML import was not building voxels because of early return if the check is disabled

[geom] Rebuild tessellated BVH after ResizeCenter

d45f1e4

[geom] Add BVH blob persistence and use it in TGeoTessellated

fcf7b4e

[geom] Support quad facets in TGeoTessellated BVH

6b18c88

Merge branch 'master' into swenzel/tessellated_on_master

30199d6

silverweed reviewed Feb 4, 2026

View reviewed changes

Add BVH-based navigation functions to TGeoTessellated #21045

Are you sure you want to change the base?

Add BVH-based navigation functions to TGeoTessellated #21045

Uh oh!

Conversation

sawenzel commented Jan 27, 2026

Uh oh!

sawenzel commented Jan 27, 2026

Uh oh!

sawenzel commented Jan 27, 2026

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

dpiparo commented Jan 28, 2026

Uh oh!

agheata left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

sawenzel commented Feb 2, 2026

Uh oh!

Uh oh!

github-actions bot commented Feb 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Test Results

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

github-actions bot commented Feb 2, 2026 •

edited

Loading