VecGeom/navigation/EmbreeNavigator.h

0001 /*
0002  * EmbreeNavigator.h
0003  *
0004  *  Created on: May 18, 2018
0005  *      Author: swenzel
0006  */
0007
0008 #ifndef VECGEOM_NAVIGATION_EMBREENAVIGATOR_H_
0009 #define VECGEOM_NAVIGATION_EMBREENAVIGATOR_H_
0010
0011 #include "VecGeom/base/Global.h"
0012
0013 #include "VecGeom/volumes/PlacedVolume.h"
0014 #include "VecGeom/base/Vector3D.h"
0015 #include "VecGeom/management/GeoManager.h"
0016 #include "VecGeom/navigation/NavigationState.h"
0017 #include "VecGeom/base/Transformation3D.h"
0018 #include "VecGeom/management/EmbreeManager.h"
0019 #include "VecGeom/navigation/VNavigator.h"
0020 #include "VecGeom/navigation/HybridSafetyEstimator.h"
0021 #include "VecGeom/navigation/SimpleABBoxNavigator.h"
0022
0023 #include <vector>
0024 #include <stack>
0025
0026 namespace vecgeom {
0027 inline namespace VECGEOM_IMPL_NAMESPACE {
0028
0029 extern double g_step;
0030 extern LogicalVolume const *g_lvol;
0031 extern VPlacedVolume const *g_pvol;
0032 extern VPlacedVolume const *g_lastexited;
0033 extern Vector3D<double> const *g_pos;
0034 extern Vector3D<double> const *g_dir;
0035 extern Vector3D<float> const *g_normals;
0036 extern int g_count;
0037 extern std::vector<int> *g_geomIDs;
0038 extern EmbreeManager::BoxIdDistancePair_t *g_hitlist;
0039
0040 // A navigator using Intel Embree as the underlying acceleration library to
0041 // exclude hit targets quickly
0042 template <bool MotherIsConvex = false>
0043 class EmbreeNavigator : public VNavigatorHelper<EmbreeNavigator<MotherIsConvex>, MotherIsConvex> {
0044
0045 private:
0046   EmbreeManager &fAccelerationManager;
0047   EmbreeNavigator()
0048       : VNavigatorHelper<EmbreeNavigator<MotherIsConvex>, MotherIsConvex>(SimpleABBoxSafetyEstimator::Instance()),
0049         fAccelerationManager(EmbreeManager::Instance())
0050   {
0051   }
0052
0053   static VPlacedVolume const *LookupDaughter(LogicalVolume const *lvol, int const daughterIndex)
0054   {
0055     return lvol->GetDaughters()[daughterIndex];
0056   }
0057
0058   // a simple sort class (based on insertionsort)
0059   template <typename T> //, typename Cmp>
0060   static void insertionsort(T *arr, unsigned int N)
0061   {
0062     for (unsigned short i = 1; i < N; ++i) {
0063       T value    = arr[i];
0064       short hole = i;
0065
0066       for (; hole > 0 && value.second < arr[hole - 1].second; --hole)
0067         arr[hole] = arr[hole - 1];
0068
0069       arr[hole] = value;
0070     }
0071   }
0072
0073   /**
0074    * Returns hitlist of daughter candidates (pairs of [daughter index, step to bounding box]) crossed by ray.
0075    */
0076   size_t GetHitCandidates_v(EmbreeManager::EmbreeAccelerationStructure const &accstructure,
0077                             Vector3D<Precision> const &point, Vector3D<Precision> const &dir,
0078                             EmbreeManager::BoxIdDistancePair_t *hitlist, float step) const
0079   {
0080     if (accstructure.fNumberObjects == 0) return false;
0081     // we need to setup an Embree ray
0082     RTCRayHit ray;
0083     ray.ray.flags = 0;
0084     ray.ray.org_x = point.x();
0085     ray.ray.org_y = point.y();
0086     ray.ray.org_z = point.z();
0087     ray.ray.dir_x = dir.x();
0088     ray.ray.dir_y = dir.y();
0089     ray.ray.dir_z = dir.z();
0090     ray.ray.tnear = 0.f;
0091     ray.ray.tfar  = 1E20f; // step is the current limit
0092
0093     g_normals = accstructure.fNormals;
0094     g_hitlist = hitlist;
0095     g_count   = 0;
0096     g_geomIDs->clear(); // = geom_seen;
0097
0098     RTCIntersectContext context;
0099     // we can't do a real capture ... so we do it via global variables
0100     auto hitFilter = [](const RTCFilterFunctionNArguments *args) {
0101       // check if a hit in this Embree structure leads to a hit
0102       // in the real geometry
0103       const auto hit = (RTCHit *)args->hit;
0104       int *hitvalid  = args->valid;
0105       const auto id  = hit->geomID;
0106       // if (g_geomIDs[id]) {
0107       //  hitvalid[0] = 0;
0108       //  return;
0109       // }
0110       // g_geomIDs[id] = true;
0111       if (std::find(g_geomIDs->begin(), g_geomIDs->end(), id) != g_geomIDs->end()) {
0112         hitvalid[0] = 0;
0113         return;
0114       }
0115       g_geomIDs->push_back(id);
0116
0117       const auto ray      = (RTCRay *)args->ray;
0118       const auto normalID = id * 6 + hit->primID;
0119       const auto normal   = g_normals[normalID];
0120       const bool backface = ray->dir_x * normal.x() + ray->dir_y * normal.y() + ray->dir_z * normal.z() > 0;
0121       float dist          = backface ? -1.f : ray->tfar;
0122
0123       // no need to sort twice !!! (in principle this thing is giving sorted inters??)
0124
0125       // std::cerr << "putting " << id << " " << dist << " " << ray->tfar << "\n";
0126       g_hitlist[g_count++] = HybridManager2::BoxIdDistancePair_t(id, dist);
0127       // we strictly take all hits
0128       hitvalid[0] = 0;
0129     };
0130
0131     rtcInitIntersectContext(&context);
0132     context.filter = hitFilter;
0133     rtcIntersect1(accstructure.fScene, &context, &ray);
0134
0135     // at this moment we have the result
0136     return g_count;
0137   }
0138
0139 public:
0140   // a generic looper function that
0141   // given an acceleration structure (an aligned bounding box hierarchy),
0142   // a hit-query will be performed, the intersected boxes sorted, looped over
0143   // and a user hook called for processing
0144   // the user hook needs to indicate with a boolean return value whether to continue looping (false)
0145   // or whether we are done (true) and can exit
0146
0147   // FIXME: might be generic enough to work for all possible kinds of BVH structures
0148   // FIXME: offer various sorting directions, etc.
0149   template <typename AccStructure, typename Func>
0150   VECGEOM_FORCE_INLINE
0151   void BVHSortedIntersectionsLooper(AccStructure const &accstructure, Vector3D<Precision> const &localpoint,
0152                                     Vector3D<Precision> const &localdir, float stepmax, Func &&userhook) const
0153   {
0154     // The following construct reserves stackspace for objects
0155     // of type IdDistPair_t WITHOUT initializing those objects
0156     using IdDistPair_t = HybridManager2::BoxIdDistancePair_t;
0157     char stackspace[VECGEOM_MAXFACETS * sizeof(IdDistPair_t)];
0158     IdDistPair_t *hitlist = reinterpret_cast<IdDistPair_t *>(&stackspace);
0159
0160     auto ncandidates = GetHitCandidates_v(accstructure, localpoint, localdir, hitlist, stepmax);
0161     // sort candidates according to their bounding volume hit distance
0162     insertionsort(hitlist, ncandidates);
0163
0164     // for (int c = 0; c < ncandidates; ++c) {
0165     //  std::cerr << "CAND " << 0 << hitlist[c].first << " " << hitlist[c].second << "\n";
0166     // }
0167
0168     for (size_t index = 0; index < ncandidates; ++index) {
0169       auto hitbox = hitlist[index];
0170       // here we got the hit candidates
0171       // now we execute user specific code to process this "hitbox"
0172       auto done = userhook(hitbox);
0173       if (done) break;
0174     }
0175   }
0176
0177   VECGEOM_FORCE_INLINE
0178   virtual bool CheckDaughterIntersections(LogicalVolume const *lvol, Vector3D<Precision> const &localpoint,
0179                                           Vector3D<Precision> const &localdir, NavigationState const *in_state,
0180                                           NavigationState * /*out_state*/, Precision &step,
0181                                           VPlacedVolume const *&hitcandidate) const override
0182   {
0183     if (lvol->GetDaughtersp()->size() == 0) return false;
0184     auto &accstructure = *fAccelerationManager.GetAccStructure(lvol);
0185
0186     BVHSortedIntersectionsLooper(
0187         accstructure, localpoint, localdir, step, [&](HybridManager2::BoxIdDistancePair_t hitbox) {
0188           // only consider those hitboxes which are within potential reach of this step
0189           if (!(step < hitbox.second)) {
0190             VPlacedVolume const *candidate = LookupDaughter(lvol, hitbox.first);
0191             Precision ddistance            = candidate->DistanceToIn(localpoint, localdir, step);
0192             const auto valid               = !IsInf(ddistance) && ddistance < step &&
0193                                !((ddistance <= 0.) && in_state && in_state->GetLastExited() == candidate);
0194             hitcandidate = valid ? candidate : hitcandidate;
0195             step         = valid ? ddistance : step;
0196             return false; // not yet done; need to continue in looper
0197           }
0198           return true; // mark done in this case
0199         });
0200     return false;
0201   }
0202
0203   //  VECGEOM_FORCE_INLINE
0204   //  virtual bool CheckDaughterIntersections(LogicalVolume const *lvol, Vector3D<Precision> const &localpoint,
0205   //                                          Vector3D<Precision> const &localdir, NavigationState const *in_state,
0206   //                                          NavigationState * /*out_state*/, Precision &step,
0207   //                                          VPlacedVolume const *&hitcandidate) const override
0208   //  {
0209   //    bool stackspace[VECGEOM_MAXFACETS * sizeof(bool)]; // alternatives are a thread_local static vector?
0210   //
0211   //    if (lvol->GetDaughtersp()->size() == 0) return false;
0212   //    const auto &accstructure = *fAccelerationManager.GetAccStructure(lvol);
0213   //
0214   //    for (int i = 0; i < lvol->GetDaughtersp()->size(); ++i) {
0215   //      stackspace[i] = false;
0216   //    }
0217   //
0218   //    // we need to setup an Embree ray
0219   //    RTCRayHit ray;
0220   //    ray.ray.flags = 0;
0221   //    ray.ray.org_x = localpoint.x();
0222   //    ray.ray.org_y = localpoint.y();
0223   //    ray.ray.org_z = localpoint.z();
0224   //    ray.ray.dir_x = localdir.x();
0225   //    ray.ray.dir_y = localdir.y();
0226   //    ray.ray.dir_z = localdir.z();
0227   //    ray.ray.tnear = 0.f;
0228   //    ray.ray.tfar  = 1E20; // step is the current limit
0229   //
0230   //    g_pos = &localpoint;
0231   //    g_dir = &localdir;
0232   //    g_step = step;
0233   //    g_lastexited = in_state ? in_state->GetLastExited() : nullptr;
0234   //    g_normals = accstructure.fNormals;
0235   //
0236   //    g_lvol = lvol;
0237   //    g_pvol = nullptr;
0238   //
0239   //    // NOT THREAD SAFE and VERY DANGEROUS; I WOULD INSTEAD VERY MUCH
0240   //    // LIKE TO BE ABLE TO CAPTURE THINGS IN THE FILTER LAMBDA
0241   //    g_geomIDs = stackspace;
0242   //    RTCIntersectContext context;
0243   //
0244   //    static int counter = 0;
0245   //    std::cerr << "# " << counter++ << "\n";
0246   //
0247   //    // we can't do a real capture ... so we do it via global variables
0248   //    auto hitFilter = [](const RTCFilterFunctionNArguments *args) {
0249   //      // check if a hit in this Embree structure leads to a hit
0250   //      // in the real geometry
0251   //
0252   //      assert(args->N == 1);
0253   //      const auto hit = (RTCHit *)args->hit;
0254   //      const auto id  = hit->geomID;
0255   //      int *hitvalid  = args->valid;
0256   //
0257   //      // if geometryID (aka bounding volume mesh) already treated we can exit
0258   //      if (g_geomIDs[id] == true) {
0259   //        hitvalid[0] = 0;
0260   //        return;
0261   //      }
0262   //      g_geomIDs[id] = true;
0263   //
0264   //      const auto ray = (RTCRay *)args->ray;
0265   //      std::cout << id << " "
0266   //                << " " << hit->primID << " " << ray->tfar << " (" << hit->Ng_x << "," << hit->Ng_y << "," <<
0267   //                hit->Ng_z
0268   //                << ") ";
0269   //      const auto normalID = id * 12 + hit->primID;
0270   //      const auto normal   = g_normals[normalID];
0271   //      const bool backface = ray->dir_x * normal.x() + ray->dir_y * normal.y() + ray->dir_z * normal.z() < 0;
0272   //      std::cout << "backface " << backface << "\n";
0273   //
0274   //      // this is assuming we get the hits in increasing distance order which somehow is not true???
0275   //      //if (!backface && g_step < ray->tfar) {
0276   //        // we can return early here (and not invalidating the hit) --> we are done!!
0277   //       // return;
0278   //      //}
0279   //
0280   //      const auto candidate = LookupDaughter(g_lvol, id);
0281   //      const auto ddistance = candidate->DistanceToIn(*g_pos, *g_dir, g_step);
0282   //      const auto valid = !IsInf(ddistance) && ddistance < g_step && !((ddistance <= 0.) && g_lastexited ==
0283   //      candidate);
0284   //      g_pvol           = valid ? candidate : g_pvol;
0285   //      g_step           = valid ? ddistance : g_step;
0286   //      hitvalid[0]      = 0;
0287   //    };
0288   //
0289   //    rtcInitIntersectContext(&context);
0290   //    context.filter = hitFilter;
0291   //    rtcIntersect1(accstructure.fScene, &context, &ray);
0292   //
0293   //    // at this moment we have the result
0294   //    hitcandidate = g_pvol;
0295   //    step = g_step;
0296   //    return false;
0297   //  }
0298
0299   static VNavigator *Instance()
0300   {
0301     static EmbreeNavigator instance;
0302     return &instance;
0303   }
0304
0305   static constexpr const char *gClassNameString = "EmbreeNavigator";
0306   typedef SimpleABBoxSafetyEstimator SafetyEstimator_t;
0307 };
0308 } // namespace VECGEOM_IMPL_NAMESPACE
0309 } // namespace vecgeom
0310
0311 #endif /* VECGEOM_NAVIGATION_EMBREENAVIGATOR_H_ */