src/INTERP_KERNEL/BBTree.txx

   1 // Copyright (C) 2007-2022  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #pragma once
  21
  22 #include <vector>
  23 #include <algorithm>
  24
  25 #include <iostream>
  26 #include <memory>
  27 #include <limits>
  28 #include <cmath>
  29
  30 constexpr double BBTREE_DFT_EPSILON = 1e-12;
  31
  32 template <int dim, class ConnType = int>
  33 class BBTree
  34 {
  35 private:
  36   BBTree* _left;
  37   BBTree* _right;
  38   int _level;
  39   double _max_left;
  40   double _min_right;
  41   const double *_bb;
  42   typename std::vector<ConnType> _elems;
  43   bool _terminal;
  44   ConnType _nbelems;
  45   double _epsilon;
  46
  47   static const int MIN_NB_ELEMS=15;
  48   static const int MAX_LEVEL=20;
  49 public:
  50   BBTree() = default;
  51   /*!
  52     Constructor of the bounding box tree
  53     \param bbs pointer to the [xmin1 xmax1 ymin1 ymax1 xmin2 xmax2 ...] array containing the bounding boxes that are to be indexed.
  54     \param elems array to the indices of the elements contained in the BBTree
  55     \param level level in the BBTree recursive structure
  56     \param nbelems nb of elements in the BBTree
  57     \param epsilon precision to which points are decided to be coincident. Epsilon can be positive or negative.
  58            If \a epsilon is positive the request method will enlarge the computed bounding box (more matching elems return).
  59            If negative the given bounding box will be tighten (less matching elems return).
  60
  61     Parameters \a elems and \a level are used only by BBTree itself for creating trees recursively. A typical use is therefore :
  62     \code
  63     int nbelems=...
  64     double* bbs= new double[2*2*nbelems];
  65     // filling bbs ...
  66     ...
  67     BBTree<2> tree = new BBTree<2>(elems,0,0,nbelems,1e-12);
  68     \endcode
  69   */
  70   BBTree(const double* bbs, ConnType* elems, int level, ConnType nbelems, double epsilon=BBTREE_DFT_EPSILON):
  71     _left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems),_epsilon(epsilon)
  72   {
  73     if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL)
  74       {
  75         _terminal=true;
  76
  77       }
  78     double median = std::numeric_limits<double>::max();
  79     {
  80       std::unique_ptr<double[]> nodes( new double [nbelems] );
  81       _elems.resize(nbelems);
  82       for (ConnType i=0; i<nbelems; i++)
  83         {
  84           ConnType elem;
  85           if (elems)
  86             elem= elems[i];
  87           else
  88             elem=i;
  89
  90           _elems[i]=elem;
  91           nodes[i]=bbs[elem*dim*2+(level%dim)*2];
  92         }
  93       if (_terminal) { return; }
  94
  95       std::nth_element<double*>(nodes.get(), nodes.get()+nbelems/2, nodes.get()+nbelems);
  96       median = *(nodes.get()+nbelems/2);
  97     }
  98     // std:: cout << *median <<std::endl;
  99
 100     std::vector<ConnType> new_elems_left;
 101     std::vector<ConnType> new_elems_right;
 102
 103     new_elems_left.reserve(nbelems/2+1);
 104     new_elems_right.reserve(nbelems/2+1);
 105     double max_left = -std::numeric_limits<double>::max();
 106     double min_right=  std::numeric_limits<double>::max();
 107     for (ConnType i=0; i<nbelems;i++)
 108       {
 109         ConnType elem;
 110         if (elems!=0)
 111           elem= elems[i];
 112         else
 113           elem=i;
 114         double max=bbs[elem*dim*2+(level%dim)*2+1];
 115         double min = bbs[elem*dim*2+(level%dim)*2];
 116
 117         if (min>median)
 118           {
 119             new_elems_right.push_back(elem);
 120             if (min<min_right) min_right = min;
 121           }
 122         else
 123
 124           {
 125             new_elems_left.push_back(elem);
 126             if (max>max_left) max_left = max;
 127           }
 128
 129
 130       }
 131     _max_left=max_left+std::abs(_epsilon);
 132     _min_right=min_right-std::abs(_epsilon);
 133     ConnType *tmp;
 134     tmp=0;
 135     if(!new_elems_left.empty())
 136       tmp=&(new_elems_left[0]);
 137     _left=new BBTree(bbs, tmp, level+1, (ConnType)new_elems_left.size(),_epsilon);
 138     tmp=0;
 139     if(!new_elems_right.empty())
 140       tmp=&(new_elems_right[0]);
 141     _right=new BBTree(bbs, tmp, level+1, (ConnType)new_elems_right.size(),_epsilon);
 142
 143   }
 144
 145   ~BBTree()
 146   {
 147     if (_left!=0)  delete _left;
 148     if (_right!=0) delete _right;
 149
 150   }
 151
 152   /*! returns in \a elems the list of elements potentially intersecting the bounding box pointed to by \a bb
 153
 154     \param bb pointer to query bounding box
 155     \param elems list of elements (given in 0-indexing that is to say in \b C \b mode) intersecting the bounding box
 156   */
 157   void getIntersectingElems(const double* bb, std::vector<ConnType>& elems) const
 158   {
 159     //  terminal node : return list of elements intersecting bb
 160     if (_terminal)
 161       {
 162         for (ConnType i=0; i<_nbelems; i++)
 163           {
 164             const double* const  bb_ptr=_bb+_elems[i]*2*dim;
 165             bool intersects = true;
 166             for (int idim=0; idim<dim; idim++)
 167               {
 168                 if (bb_ptr[idim*2]-bb[idim*2+1]>-_epsilon|| bb_ptr[idim*2+1]-bb[idim*2]<_epsilon)
 169                   intersects=false;
 170               }
 171             if (intersects)
 172               {
 173                 elems.push_back(_elems[i]);
 174               }
 175           }
 176         return;
 177       }
 178
 179     //non terminal node
 180     double min = bb[(_level%dim)*2];
 181     double max = bb[(_level%dim)*2+1];
 182     if (max < _min_right)
 183       {
 184         _left->getIntersectingElems(bb, elems);
 185         return;
 186       }
 187     if (min> _max_left)
 188       {
 189         _right->getIntersectingElems(bb,elems);
 190         return;
 191       }
 192     _left->getIntersectingElems(bb,elems);
 193     _right->getIntersectingElems(bb,elems);
 194   }
 195
 196   /*!
 197    * This method is very close to getIntersectingElems except that it returns number of elems instead of elems themselves.
 198    */
 199   ConnType getNbOfIntersectingElems(const double* bb)
 200   {
 201     //  terminal node : return list of elements intersecting bb
 202     ConnType ret(0);
 203     if (_terminal)
 204       {
 205         for (ConnType i=0; i<_nbelems; i++)
 206           {
 207             const double* const  bb_ptr=_bb+_elems[i]*2*dim;
 208             bool intersects = true;
 209             for (int idim=0; idim<dim; idim++)
 210               {
 211                 if (bb_ptr[idim*2]-bb[idim*2+1]>-_epsilon|| bb_ptr[idim*2+1]-bb[idim*2]<_epsilon)
 212                   intersects=false;
 213               }
 214             if (intersects)
 215               ret++;
 216           }
 217         return ret;
 218       }
 219     //non terminal node
 220     double min = bb[(_level%dim)*2];
 221     double max = bb[(_level%dim)*2+1];
 222     if (max < _min_right)
 223       return _left->getNbOfIntersectingElems(bb);
 224     if (min> _max_left)
 225       return _right->getNbOfIntersectingElems(bb);
 226     return _left->getNbOfIntersectingElems(bb)+_right->getNbOfIntersectingElems(bb);
 227   }
 228
 229
 230   /*! returns in \a elems the list of elements potentially containing the point pointed to by \a xx
 231     \param xx pointer to query point coords
 232     \param elems list of elements (given in 0-indexing) intersecting the bounding box
 233   */
 234   void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const
 235   {
 236     //  terminal node : return list of elements intersecting bb
 237     if (_terminal)
 238       {
 239         for (ConnType i=0; i<_nbelems; i++)
 240           {
 241             const double* const  bb_ptr=_bb+_elems[i]*2*dim;
 242             bool intersects = true;
 243             for (int idim=0; idim<dim; idim++)
 244               {
 245                 if (bb_ptr[idim*2]-xx[idim]>_epsilon|| bb_ptr[idim*2+1]-xx[idim]<-_epsilon)
 246                   intersects=false;
 247               }
 248             if (intersects)
 249               {
 250                 elems.push_back(_elems[i]);
 251               }
 252           }
 253         return;
 254       }
 255
 256     //non terminal node
 257     if (xx[_level%dim] < _min_right)
 258       {
 259         _left->getElementsAroundPoint(xx, elems);
 260         return;
 261       }
 262     if (xx[_level%dim]> _max_left)
 263       {
 264         _right->getElementsAroundPoint(xx,elems);
 265         return;
 266       }
 267     _left->getElementsAroundPoint(xx,elems);
 268     _right->getElementsAroundPoint(xx,elems);
 269   }
 270
 271   ConnType size()
 272   {
 273     if (_terminal) return _nbelems;
 274     return _left->size()+_right->size();
 275   }
 276 };