include/Geant4/G4KDTree.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // Author: Mathieu Karamitros
0028
0029 // The code is developed in the framework of the ESA AO7146
0030 //
0031 // We would be very happy hearing from you, send us your feedback! :)
0032 //
0033 // In order for Geant4-DNA to be maintained and still open-source,
0034 // article citations are crucial.
0035 // If you use Geant4-DNA chemistry and you publish papers about your software,
0036 // in addition to the general paper on Geant4-DNA:
0037 //
0038 // Int. J. Model. Simul. Sci. Comput. 1 (2010) 157–178
0039 //
0040 // we would be very happy if you could please also cite the following
0041 // reference papers on chemistry:
0042 //
0043 // J. Comput. Phys. 274 (2014) 841-882
0044 // Prog. Nucl. Sci. Tec. 2 (2011) 503-508
0045
0046 #ifndef G4KDTREE_HH
0047 #define G4KDTREE_HH 1
0048
0049 #include <vector>
0050 #include "G4Types.hh"
0051 #include "G4KDNode.hh"
0052 #include "G4KDTreeResult.hh"
0053
0054 class G4KDMap;
0055 template <typename PointT>
0056 class G4KDNode;
0057
0058 //__________________________________
0059 // Methods to act on kdnode
0060 // Methods defined in G4KDNode.cc :
0061 void InactiveNode(G4KDNode_Base*);
0062 void Free(G4KDNode_Base*&);
0063 //__________________________________
0064
0065 /**
0066  * G4KDTree is used by the ITManager to locate the neareast neighbours.
0067  * A kdtree sorts out node in such a way that it reduces the number of node
0068  * check. The results of this search can be retrieved by G4KDTreeResultHandle.
0069  */
0070 class G4KDTree
0071 {
0072   friend class G4KDNode_Base;
0073
0074  public:
0075   G4KDTree(std::size_t dim = 3);
0076   ~G4KDTree();
0077   void Clear();
0078
0079   void Print(std::ostream& out = G4cout) const;
0080   void Build();
0081   void NoticeNodeDeactivation()
0082   {
0083     fNbActiveNodes--;
0084     if(fNbActiveNodes <= 0)
0085     {
0086       Clear();
0087     }
0088   }
0089
0090   std::size_t GetDim() const { return fDim; }
0091   G4int GetNbNodes() const { return fNbNodes; }
0092   G4KDNode_Base* GetRoot() { return fRoot; }
0093
0094   template <typename PointT>
0095   G4KDNode_Base* InsertMap(PointT* pos);
0096
0097   // Insert and attache the data to a node at the specified position
0098   // In return, it gives you the corresponding node
0099   template <typename PointT>
0100   G4KDNode_Base* Insert(PointT* pos);  // 3D
0101
0102   template <typename PointT>
0103   G4KDNode_Base* Insert(const PointT& pos);  // 3D
0104
0105   /* Find one of the nearest nodes from the specified point.
0106    *
0107    * This function returns a pointer to a result set with at most one element.
0108    */
0109   template <typename Position>
0110   G4KDTreeResultHandle Nearest(const Position& pos);
0111   G4KDTreeResultHandle Nearest(G4KDNode_Base* node);
0112
0113   /* Find any nearest nodes from the specified point within a range.
0114    *
0115    * This function returns a pointer to a result set, which can be manipulated
0116    * by the G4KDTreeResult.
0117    * The returned pointer can be null as an indication of an error. Otherwise
0118    * a valid result set is always returned which may contain 0 or more elements.
0119    */
0120   template <typename Position>
0121   G4KDTreeResultHandle NearestInRange(const Position& pos, const G4double& range);
0122   G4KDTreeResultHandle NearestInRange(G4KDNode_Base* node, const G4double& range);
0123
0124   void* operator new(std::size_t);
0125   void operator delete(void*);
0126
0127  protected:
0128   //______________________________________________________________________
0129   class HyperRect
0130   {
0131    public:
0132     HyperRect(std::size_t dim)
0133       : fDim(dim)
0134       , fMin(new G4double[fDim])
0135       , fMax(new G4double[fDim])
0136     {}
0137
0138     template <typename Position>
0139     void SetMinMax(const Position& min, const Position& max)
0140     {
0141       for(std::size_t i = 0; i < fDim; ++i)
0142       {
0143         fMin[i] = min[(G4int)i];
0144         fMax[i] = max[(G4int)i];
0145       }
0146     }
0147
0148     ~HyperRect()
0149     {
0150       delete[] fMin;
0151       delete[] fMax;
0152     }
0153
0154     HyperRect(const HyperRect& rect)
0155     {
0156       fDim = rect.fDim;
0157       fMin = new G4double[fDim];
0158       fMax = new G4double[fDim];
0159
0160       for(std::size_t i = 0; i < fDim; ++i)
0161       {
0162         fMin[i] = rect.fMin[i];
0163         fMax[i] = rect.fMax[i];
0164       }
0165     }
0166
0167     template <typename Position>
0168     void Extend(const Position& pos)
0169     {
0170       for(G4int i = 0; i < (G4int)fDim; ++i)
0171       {
0172         if(pos[i] < fMin[i])
0173         {
0174           fMin[i] = pos[i];
0175         }
0176         if(pos[i] > fMax[i])
0177         {
0178           fMax[i] = pos[i];
0179         }
0180       }
0181     }
0182
0183     template <typename Position>
0184     G4bool CompareDistSqr(const Position& pos, const G4double* bestmatch)
0185     {
0186       G4double result = 0;
0187
0188       for(std::size_t i = 0; i < fDim; ++i)
0189       {
0190         if(pos[(G4int)i] < fMin[i])
0191         {
0192           result += sqr(fMin[i] - pos[(G4int)i]);
0193         }
0194         else if(pos[(G4int)i] > fMax[i])
0195         {
0196           result += sqr(fMax[i] - pos[(G4int)i]);
0197         }
0198
0199         if(result >= *bestmatch){
0200           return false;
0201         }
0202       }
0203
0204       return true;
0205     }
0206
0207     std::size_t GetDim() { return fDim; }
0208     G4double* GetMin() { return fMin; }
0209     G4double* GetMax() { return fMax; }
0210
0211    protected:
0212     std::size_t fDim;
0213     G4double *fMin, *fMax; /* minimum/maximum coords */
0214
0215    private:
0216     // should not be used
0217     HyperRect& operator=(const HyperRect& rhs)
0218     {
0219       if (this == &rhs) return *this;
0220       return *this;
0221     }
0222   };
0223
0224  protected:
0225   void __InsertMap(G4KDNode_Base* node);
0226   void __Clear_Rec(G4KDNode_Base* node);
0227
0228   template <typename Position>
0229   G4int __NearestInRange(G4KDNode_Base* node, const Position& pos,
0230                        const G4double& range_sq, const G4double& range,
0231                        G4KDTreeResult& list, G4int ordered,
0232                        G4KDNode_Base* source_node = nullptr);
0233
0234   template <typename Position>
0235   void __NearestToPosition(G4KDNode_Base* node, const Position& pos,
0236                            G4KDNode_Base*& result, G4double* result_dist_sq,
0237                            HyperRect* fRect);
0238
0239   template <typename Position>
0240   void __NearestToNode(G4KDNode_Base* source_node, G4KDNode_Base* node,
0241                        const Position& pos, std::vector<G4KDNode_Base*>& result,
0242                        G4double* result_dist_sq, HyperRect* fRect, G4int& nbresult);
0243
0244  protected:
0245   HyperRect* fRect     = nullptr;
0246   G4KDNode_Base* fRoot = nullptr;
0247   std::size_t fDim;
0248   G4int fNbNodes       = 0;
0249   G4int fNbActiveNodes = 0;
0250   G4KDMap* fKDMap;
0251   static G4Allocator<G4KDTree>*& fgAllocator();
0252 };
0253
0254 #include "G4KDTree.icc"
0255
0256 #endif  // G4KDTREE_HH