src/INTERP_KERNEL/BBTree.txx

   1 // Copyright (C) 2007-2014  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 #ifndef __BBTREE_TXX__
  20 #define __BBTREE_TXX__
  21
  22 #include <vector>
  23 #include <algorithm>
  24
  25 #include <iostream>
  26 #include <limits>
  27 #include <cmath>
  28
  29 template <int dim, class ConnType = int>
  30 class BBTree
  31 {
  32
  33 private:
  34   BBTree* _left;
  35   BBTree* _right;
  36   int _level;
  37   double _max_left;
  38   double _min_right;
  39   const double *_bb;
  40   typename std::vector<ConnType> _elems;
  41   bool _terminal;
  42   ConnType _nbelems;
  43   double _epsilon;
  44
  45   static const int MIN_NB_ELEMS=15;
  46   static const int MAX_LEVEL=20;
  47 public:
  48
  49   /*!
  50     Constructor of the bounding box tree
  51     \param bbs pointer to the [xmin1 xmax1 ymin1 ymax1 xmin2 xmax2 ...] array containing the bounding boxes that are to be indexed.
  52     \param elems array to the indices of the elements contained in the BBTree
  53     \param level level in the BBTree recursive structure
  54     \param nbelems nb of elements in the BBTree
  55     \param epsilon precision to which points are decided to be coincident. Epsilon can be positive or negative.
  56            If \a epsilon is positive the request method will enlarge the computed bounding box (more matching elems return).
  57            If negative the given bounding box will be tighten (less matching elems return).
  58
  59     Parameters \a elems and \a level are used only by BBTree itself for creating trees recursively. A typical use is therefore :
  60     \code
  61     int nbelems=...
  62     double* bbs= new double[2*2*nbelems];
  63     // filling bbs ...
  64     ...
  65     BBTree<2> tree = new BBTree<2>(elems,0,0,nbelems,1e-12);
  66     \endcode
  67   */
  68
  69   BBTree(const double* bbs, ConnType* elems, int level, ConnType nbelems, double epsilon=1e-12):
  70     _left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems),_epsilon(epsilon)
  71   {
  72     if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL)
  73       {
  74         _terminal=true;
  75
  76       }
  77     double* nodes=new double [nbelems];
  78     _elems.resize(nbelems);
  79     for (ConnType i=0; i<nbelems; i++)
  80       {
  81         ConnType elem;
  82         if (elems!=0)
  83           elem= elems[i];
  84         else
  85           elem=i;
  86
  87         _elems[i]=elem;
  88         nodes[i]=bbs[elem*dim*2+(level%dim)*2];
  89       }
  90     if (_terminal) { delete[] nodes; return;}
  91
  92     std::nth_element<double*>(nodes, nodes+nbelems/2, nodes+nbelems);
  93     double median = *(nodes+nbelems/2);
  94     delete[] nodes;
  95     // std:: cout << *median <<std::endl;
  96
  97     std::vector<ConnType> new_elems_left;
  98     std::vector<ConnType> new_elems_right;
  99
 100     new_elems_left.reserve(nbelems/2+1);
 101     new_elems_right.reserve(nbelems/2+1);
 102     double max_left = -std::numeric_limits<double>::max();
 103     double min_right=  std::numeric_limits<double>::max();
 104     for (int i=0; i<nbelems;i++)
 105       {
 106         int elem;
 107         if (elems!=0)
 108           elem= elems[i];
 109         else
 110           elem=i;
 111         double max=bbs[elem*dim*2+(level%dim)*2+1];
 112         double min = bbs[elem*dim*2+(level%dim)*2];
 113
 114         if (min>median)
 115           {
 116             new_elems_right.push_back(elem);
 117             if (min<min_right) min_right = min;
 118           }
 119         else
 120
 121           {
 122             new_elems_left.push_back(elem);
 123             if (max>max_left) max_left = max;
 124           }
 125
 126
 127       }
 128     _max_left=max_left+std::abs(_epsilon);
 129     _min_right=min_right-std::abs(_epsilon);
 130     ConnType *tmp;
 131     tmp=0;
 132     if(!new_elems_left.empty())
 133       tmp=&(new_elems_left[0]);
 134     _left=new BBTree(bbs, tmp, level+1, (int)new_elems_left.size(),_epsilon);
 135     tmp=0;
 136     if(!new_elems_right.empty())
 137       tmp=&(new_elems_right[0]);
 138     _right=new BBTree(bbs, tmp, level+1, (int)new_elems_right.size(),_epsilon);
 139
 140   }
 141
 142
 143
 144   ~BBTree()
 145   {
 146     if (_left!=0)  delete _left;
 147     if (_right!=0) delete _right;
 148
 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
 158   void getIntersectingElems(const double* bb, std::vector<ConnType>& elems) const
 159   {
 160     //  terminal node : return list of elements intersecting bb
 161     if (_terminal)
 162       {
 163         for (int i=0; i<_nbelems; i++)
 164           {
 165             const double* const  bb_ptr=_bb+_elems[i]*2*dim;
 166             bool intersects = true;
 167             for (int idim=0; idim<dim; idim++)
 168               {
 169                 if (bb_ptr[idim*2]-bb[idim*2+1]>-_epsilon|| bb_ptr[idim*2+1]-bb[idim*2]<_epsilon)
 170                   intersects=false;
 171               }
 172             if (intersects)
 173               {
 174                 elems.push_back(_elems[i]);
 175               }
 176           }
 177         return;
 178       }
 179
 180     //non terminal node
 181     double min = bb[(_level%dim)*2];
 182     double max = bb[(_level%dim)*2+1];
 183     if (max < _min_right)
 184       {
 185         _left->getIntersectingElems(bb, elems);
 186         return;
 187       }
 188     if (min> _max_left)
 189       {
 190         _right->getIntersectingElems(bb,elems);
 191         return;
 192       }
 193     _left->getIntersectingElems(bb,elems);
 194     _right->getIntersectingElems(bb,elems);
 195   }
 196
 197   /*!
 198    * This method is very close to getIntersectingElems except that it returns number of elems instead of elems themselves.
 199    */
 200   int getNbOfIntersectingElems(const double* bb)
 201   {
 202     //  terminal node : return list of elements intersecting bb
 203     int ret(0);
 204     if (_terminal)
 205       {
 206         for (int i=0; i<_nbelems; i++)
 207           {
 208             const double* const  bb_ptr=_bb+_elems[i]*2*dim;
 209             bool intersects = true;
 210             for (int idim=0; idim<dim; idim++)
 211               {
 212                 if (bb_ptr[idim*2]-bb[idim*2+1]>-_epsilon|| bb_ptr[idim*2+1]-bb[idim*2]<_epsilon)
 213                   intersects=false;
 214               }
 215             if (intersects)
 216               ret++;
 217           }
 218         return ret;
 219       }
 220     //non terminal node
 221     double min = bb[(_level%dim)*2];
 222     double max = bb[(_level%dim)*2+1];
 223     if (max < _min_right)
 224       return _left->getNbOfIntersectingElems(bb);
 225     if (min> _max_left)
 226       return _right->getNbOfIntersectingElems(bb);
 227     return _left->getNbOfIntersectingElems(bb)+_right->getNbOfIntersectingElems(bb);
 228   }
 229
 230
 231   /*! returns in \a elems the list of elements potentially containing the point pointed to by \a xx
 232     \param xx pointer to query point coords
 233     \param elems list of elements (given in 0-indexing) intersecting the bounding box
 234   */
 235
 236   void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const
 237   {
 238     //  terminal node : return list of elements intersecting bb
 239     if (_terminal)
 240       {
 241         for (ConnType i=0; i<_nbelems; i++)
 242           {
 243             const double* const  bb_ptr=_bb+_elems[i]*2*dim;
 244             bool intersects = true;
 245             for (int idim=0; idim<dim; idim++)
 246               {
 247                 if (bb_ptr[idim*2]-xx[idim]>_epsilon|| bb_ptr[idim*2+1]-xx[idim]<-_epsilon)
 248                   intersects=false;
 249               }
 250             if (intersects)
 251               {
 252                 elems.push_back(_elems[i]);
 253               }
 254           }
 255         return;
 256       }
 257
 258     //non terminal node
 259     if (xx[_level%dim] < _min_right)
 260       {
 261         _left->getElementsAroundPoint(xx, elems);
 262         return;
 263       }
 264     if (xx[_level%dim]> _max_left)
 265       {
 266         _right->getElementsAroundPoint(xx,elems);
 267         return;
 268       }
 269     _left->getElementsAroundPoint(xx,elems);
 270     _right->getElementsAroundPoint(xx,elems);
 271   }
 272
 273
 274
 275   int size()
 276   {
 277     if (_terminal) return _nbelems;
 278     return _left->size()+_right->size();
 279   }
 280 };
 281 #endif