VecGeom/management/HybridManager2.h

0001 /*
0002  * HybridManager2.h
0003  *
0004  *  Created on: 27.08.2015 by yang.zhang@cern.ch
0005  *  integrated into main development line by sandro.wenzel@cern.ch 24.11.2015
0006  *
0007  */
0008
0009 #ifndef VECGEOM_HYBRIDMANAGER_H
0010 #define VECGEOM_HYBRIDMANAGER_H
0011
0012 #include "VecGeom/base/Global.h"
0013
0014 #include "VecGeom/volumes/PlacedVolume.h"
0015 #include "VecGeom/base/Vector3D.h"
0016 #include "VecGeom/management/GeoManager.h"
0017 #include "VecGeom/base/Transformation3D.h"
0018 #include "VecGeom/volumes/kernel/BoxImplementation.h"
0019 #include "VecGeom/base/AlignmentAllocator.h"
0020 #include "VecGeom/management/ABBoxManager.h"
0021
0022 #include <queue>
0023 #include <map>
0024 #include <vector>
0025
0026 namespace vecgeom {
0027 inline namespace VECGEOM_IMPL_NAMESPACE {
0028
0029 // a singleton class which manages a particular helper structure for voxelized navigation
0030 // the helper structure is a hybrid between a flat list of aligned bounding boxed and a pure boundary volume hierarchy
0031 // (BVH)
0032 // and is hence called "HybridManager": TODO: come up with a more appropriate name
0033 class HybridManager2 {
0034
0035 public:
0036   using Float_v = vecgeom::VectorBackend::Float_v;
0037   typedef float Real_t;
0038   typedef Vector3D<Float_v> ABBox_v;
0039
0040   // scalar or vector vectors
0041   typedef Vector3D<Precision> ABBox_s;
0042   // use old style arrays here as std::vector has some problems
0043   // with Vector3D<kVc::Double_t>
0044   typedef ABBox_s *ABBoxContainer_t;
0045   typedef ABBox_v *ABBoxContainer_v;
0046
0047   // first index is # daughter index, second is step
0048   typedef std::pair<int, double> BoxIdDistancePair_t;
0049   using HitContainer_t = std::vector<BoxIdDistancePair_t>;
0050
0051   // the actual class encapsulating the bounding boxes + tree structure information
0052   // for a shallow bounding volume hierarchy acceleration structure
0053   struct HybridBoxAccelerationStructure {
0054     size_t fNumberOfOriginalBoxes = 0; // the number of objects/original bounding boxes this structure is representing
0055     ABBoxContainer_v fABBoxes_v   = nullptr;
0056     std::vector<int> *fNodeToDaughters = nullptr;
0057   };
0058
0059 private:
0060   // keeps/registers an acceleration structure for logical volumes
0061   std::vector<HybridBoxAccelerationStructure const *> fStructureHolder;
0062
0063 public:
0064   // initialized the helper structure for a given logical volume
0065   void InitStructure(LogicalVolume const *lvol);
0066
0067   // initialized the helper structure for the complete geometry
0068   void InitVoxelStructureForCompleteGeometry()
0069   {
0070     std::vector<LogicalVolume const *> logicalvolumes;
0071     GeoManager::Instance().GetAllLogicalVolumes(logicalvolumes);
0072     for (auto lvol : logicalvolumes) {
0073       InitStructure(lvol);
0074     }
0075   }
0076
0077   static HybridManager2 &Instance()
0078   {
0079     static HybridManager2 manager;
0080     return manager;
0081   }
0082
0083   // removed/deletes the helper structure for a given logical volume
0084   void RemoveStructure(LogicalVolume const *lvol);
0085
0086   template <typename C, typename Compare>
0087   static void sort(C &v, Compare cmp)
0088   {
0089     std::sort(v.begin(), v.end(), cmp);
0090   }
0091
0092   // verbose output of helper structure
0093   VPlacedVolume const *PrintHybrid(LogicalVolume const *) const;
0094
0095   ABBoxContainer_v GetABBoxes_v(HybridBoxAccelerationStructure const &structure, int &size, int &numberOfNodes) const
0096   {
0097     constexpr auto kVS = vecCore::VectorSize<Float_v>();
0098     assert(structure.fNumberOfOriginalBoxes != 0);
0099     int numberOfFirstLevelNodes =
0100         structure.fNumberOfOriginalBoxes / kVS + (structure.fNumberOfOriginalBoxes % kVS == 0 ? 0 : 1);
0101     numberOfNodes = numberOfFirstLevelNodes + structure.fNumberOfOriginalBoxes;
0102     size = numberOfFirstLevelNodes / kVS + (numberOfFirstLevelNodes % kVS == 0 ? 0 : 1) + numberOfFirstLevelNodes;
0103     assert(structure.fABBoxes_v != nullptr);
0104     return structure.fABBoxes_v;
0105   }
0106
0107   HybridBoxAccelerationStructure const *GetAccStructure(LogicalVolume const *lvol) const
0108   {
0109     return fStructureHolder[lvol->id()];
0110   }
0111
0112   // public method allowing to build a hybrid acceleration structure
0113   // given a vector of aligned bounding boxes
0114   // can be used by clients (not coupled to LogicalVolumes)
0115   HybridBoxAccelerationStructure *BuildStructure(ABBoxManager::ABBoxContainer_t alignedboxes,
0116                                                  size_t numberofboxes) const;
0117
0118 private:
0119   // private methods use
0120
0121   template <typename Container_t>
0122   void EqualizeClusters(Container_t &clusters, SOA3D<Precision> &centers, SOA3D<Precision> const &allvolumecenters,
0123                         size_t const maxNodeSize);
0124
0125   template <typename Container_t>
0126   void InitClustersWithKMeans(ABBoxManager::ABBoxContainer_t, int, Container_t &, SOA3D<Precision> &,
0127                               SOA3D<Precision> &, int const numberOfInterations = 50) const;
0128
0129   void RecalculateCentres(SOA3D<Precision> &centers, SOA3D<Precision> const &allvolumecenters,
0130                           std::vector<std::vector<int>> const &clusters);
0131   struct OrderByDistance {
0132     bool operator()(std::pair<int, Precision> const &a, std::pair<int, Precision> const &b)
0133     {
0134       return a.second < b.second;
0135     }
0136   };
0137   typedef std::priority_queue<std::pair<int, Precision>, std::vector<std::pair<int, Precision>>, OrderByDistance>
0138       distanceQueue;
0139   void AssignVolumesToClusters(std::vector<std::vector<int>> &clusters, SOA3D<Precision> const &centers,
0140                                SOA3D<Precision> const &allvolumecenters);
0141
0142   void BuildStructure_v(LogicalVolume const *vol);
0143
0144   static bool IsBiggerCluster(std::vector<int> const &first, std::vector<int> const &second)
0145   {
0146     return first.size() > second.size();
0147   }
0148
0149 }; // end class
0150 }
0151 } // end namespace
0152
0153 #endif