src/INTERP_KERNEL/BBTree.txx

   1 // Copyright (C) 2007-2019  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 <limits>
  27 #include <cmath>
  28
  29 constexpr double BBTREE_DFT_EPSILON = 1e-12;
  30
  31 template <int dim, class ConnType = int>
  32 class BBTree
  33 {
  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
  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
  71   BBTree(const double* bbs, ConnType* elems, int level, ConnType nbelems, double epsilon=BBTREE_DFT_EPSILON):
  72     _left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems),_epsilon(epsilon)
  73   {
  74     if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL)
  75       {
  76         _terminal=true;
  77
  78       }
  79     double* nodes=new double [nbelems];
  80     _elems.resize(nbelems);
  81     for (ConnType i=0; i<nbelems; i++)
  82       {
  83         ConnType elem;
  84         if (elems!=0)
  85           elem= elems[i];
  86         else
  87           elem=i;
  88
  89         _elems[i]=elem;
  90         nodes[i]=bbs[elem*dim*2+(level%dim)*2];
  91       }
  92     if (_terminal) { delete[] nodes; return;}
  93
  94     std::nth_element<double*>(nodes, nodes+nbelems/2, nodes+nbelems);
  95     double median = *(nodes+nbelems/2);
  96     delete[] nodes;
  97     // std:: cout << *median <<std::endl;
  98
  99     std::vector<ConnType> new_elems_left;
 100     std::vector<ConnType> new_elems_right;
 101
 102     new_elems_left.reserve(nbelems/2+1);
 103     new_elems_right.reserve(nbelems/2+1);
 104     double max_left = -std::numeric_limits<double>::max();
 105     double min_right=  std::numeric_limits<double>::max();
 106     for (ConnType i=0; i<nbelems;i++)
 107       {
 108         ConnType elem;
 109         if (elems!=0)
 110           elem= elems[i];
 111         else
 112           elem=i;
 113         double max=bbs[elem*dim*2+(level%dim)*2+1];
 114         double min = bbs[elem*dim*2+(level%dim)*2];
 115
 116         if (min>median)
 117           {
 118             new_elems_right.push_back(elem);
 119             if (min<min_right) min_right = min;
 120           }
 121         else
 122
 123           {
 124             new_elems_left.push_back(elem);
 125             if (max>max_left) max_left = max;
 126           }
 127
 128
 129       }
 130     _max_left=max_left+std::abs(_epsilon);
 131     _min_right=min_right-std::abs(_epsilon);
 132     ConnType *tmp;
 133     tmp=0;
 134     if(!new_elems_left.empty())
 135       tmp=&(new_elems_left[0]);
 136     _left=new BBTree(bbs, tmp, level+1, (ConnType)new_elems_left.size(),_epsilon);
 137     tmp=0;
 138     if(!new_elems_right.empty())
 139       tmp=&(new_elems_right[0]);
 140     _right=new BBTree(bbs, tmp, level+1, (ConnType)new_elems_right.size(),_epsilon);
 141
 142   }
 143
 144
 145
 146   ~BBTree()
 147   {
 148     if (_left!=0)  delete _left;
 149     if (_right!=0) delete _right;
 150
 151   }
 152
 153
 154   /*! returns in \a elems the list of elements potentially intersecting the bounding box pointed to by \a bb
 155
 156     \param bb pointer to query bounding box
 157     \param elems list of elements (given in 0-indexing that is to say in \b C \b mode) intersecting the bounding box
 158   */
 159
 160   void getIntersectingElems(const double* bb, std::vector<ConnType>& elems) const
 161   {
 162     //  terminal node : return list of elements intersecting bb
 163     if (_terminal)
 164       {
 165         for (ConnType i=0; i<_nbelems; i++)
 166           {
 167             const double* const  bb_ptr=_bb+_elems[i]*2*dim;
 168             bool intersects = true;
 169             for (int idim=0; idim<dim; idim++)
 170               {
 171                 if (bb_ptr[idim*2]-bb[idim*2+1]>-_epsilon|| bb_ptr[idim*2+1]-bb[idim*2]<_epsilon)
 172                   intersects=false;
 173               }
 174             if (intersects)
 175               {
 176                 elems.push_back(_elems[i]);
 177               }
 178           }
 179         return;
 180       }
 181
 182     //non terminal node
 183     double min = bb[(_level%dim)*2];
 184     double max = bb[(_level%dim)*2+1];
 185     if (max < _min_right)
 186       {
 187         _left->getIntersectingElems(bb, elems);
 188         return;
 189       }
 190     if (min> _max_left)
 191       {
 192         _right->getIntersectingElems(bb,elems);
 193         return;
 194       }
 195     _left->getIntersectingElems(bb,elems);
 196     _right->getIntersectingElems(bb,elems);
 197   }
 198
 199   /*!
 200    * This method is very close to getIntersectingElems except that it returns number of elems instead of elems themselves.
 201    */
 202   ConnType getNbOfIntersectingElems(const double* bb)
 203   {
 204     //  terminal node : return list of elements intersecting bb
 205     ConnType ret(0);
 206     if (_terminal)
 207       {
 208         for (ConnType i=0; i<_nbelems; i++)
 209           {
 210             const double* const  bb_ptr=_bb+_elems[i]*2*dim;
 211             bool intersects = true;
 212             for (int idim=0; idim<dim; idim++)
 213               {
 214                 if (bb_ptr[idim*2]-bb[idim*2+1]>-_epsilon|| bb_ptr[idim*2+1]-bb[idim*2]<_epsilon)
 215                   intersects=false;
 216               }
 217             if (intersects)
 218               ret++;
 219           }
 220         return ret;
 221       }
 222     //non terminal node
 223     double min = bb[(_level%dim)*2];
 224     double max = bb[(_level%dim)*2+1];
 225     if (max < _min_right)
 226       return _left->getNbOfIntersectingElems(bb);
 227     if (min> _max_left)
 228       return _right->getNbOfIntersectingElems(bb);
 229     return _left->getNbOfIntersectingElems(bb)+_right->getNbOfIntersectingElems(bb);
 230   }
 231
 232
 233   /*! returns in \a elems the list of elements potentially containing the point pointed to by \a xx
 234     \param xx pointer to query point coords
 235     \param elems list of elements (given in 0-indexing) intersecting the bounding box
 236   */
 237
 238   void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const
 239   {
 240     //  terminal node : return list of elements intersecting bb
 241     if (_terminal)
 242       {
 243         for (ConnType i=0; i<_nbelems; i++)
 244           {
 245             const double* const  bb_ptr=_bb+_elems[i]*2*dim;
 246             bool intersects = true;
 247             for (int idim=0; idim<dim; idim++)
 248               {
 249                 if (bb_ptr[idim*2]-xx[idim]>_epsilon|| bb_ptr[idim*2+1]-xx[idim]<-_epsilon)
 250                   intersects=false;
 251               }
 252             if (intersects)
 253               {
 254                 elems.push_back(_elems[i]);
 255               }
 256           }
 257         return;
 258       }
 259
 260     //non terminal node
 261     if (xx[_level%dim] < _min_right)
 262       {
 263         _left->getElementsAroundPoint(xx, elems);
 264         return;
 265       }
 266     if (xx[_level%dim]> _max_left)
 267       {
 268         _right->getElementsAroundPoint(xx,elems);
 269         return;
 270       }
 271     _left->getElementsAroundPoint(xx,elems);
 272     _right->getElementsAroundPoint(xx,elems);
 273   }
 274
 275
 276
 277   ConnType size()
 278   {
 279     if (_terminal) return _nbelems;
 280     return _left->size()+_right->size();
 281   }
 282 };