src/INTERP_KERNEL/BBTree.txx

   1 // Copyright (C) 2007-2012  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.
   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
  56
  57     Parameters \a elems and \a level are used only by BBTree itself for creating trees recursively. A typical use is therefore :
  58     \code
  59     int nbelems=...
  60     double* bbs= new double[2*2*nbelems];
  61     // filling bbs ...
  62     ...
  63     BBTree<2> tree = new BBTree<2>(elems,0,0,nbelems,1e-12);
  64     \endcode
  65   */
  66
  67   BBTree(const double* bbs, ConnType* elems, int level, ConnType nbelems, double epsilon=1E-12):
  68     _left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems),_epsilon(epsilon)
  69   {
  70     if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL)
  71       {
  72         _terminal=true;
  73
  74       }
  75     double* nodes=new double [nbelems];
  76     _elems.resize(nbelems);
  77     for (ConnType i=0; i<nbelems; i++)
  78       {
  79         ConnType elem;
  80         if (elems!=0)
  81           elem= elems[i];
  82         else
  83           elem=i;
  84
  85         _elems[i]=elem;
  86         nodes[i]=bbs[elem*dim*2+(level%dim)*2];
  87       }
  88     if (_terminal) { delete[] nodes; return;}
  89
  90     std::nth_element<double*>(nodes, nodes+nbelems/2, nodes+nbelems);
  91     double median = *(nodes+nbelems/2);
  92     delete[] nodes;
  93     // std:: cout << *median <<std::endl;
  94
  95     std::vector<ConnType> new_elems_left;
  96     std::vector<ConnType> new_elems_right;
  97
  98     new_elems_left.reserve(nbelems/2+1);
  99     new_elems_right.reserve(nbelems/2+1);
 100     double max_left = -std::numeric_limits<double>::max();
 101     double min_right=  std::numeric_limits<double>::max();
 102     for (int i=0; i<nbelems;i++)
 103       {
 104         int elem;
 105         if (elems!=0)
 106           elem= elems[i];
 107         else
 108           elem=i;
 109         double max=bbs[elem*dim*2+(level%dim)*2+1];
 110         double min = bbs[elem*dim*2+(level%dim)*2];
 111
 112         if (min>median)
 113           {
 114             new_elems_right.push_back(elem);
 115             if (min<min_right) min_right = min;
 116           }
 117         else
 118
 119           {
 120             new_elems_left.push_back(elem);
 121             if (max>max_left) max_left = max;
 122           }
 123
 124
 125       }
 126     _max_left=max_left+std::abs(_epsilon);
 127     _min_right=min_right-std::abs(_epsilon);
 128     ConnType *tmp;
 129     tmp=0;
 130     if(!new_elems_left.empty())
 131       tmp=&(new_elems_left[0]);
 132     _left=new BBTree(bbs, tmp, level+1, (int)new_elems_left.size(),_epsilon);
 133     tmp=0;
 134     if(!new_elems_right.empty())
 135       tmp=&(new_elems_right[0]);
 136     _right=new BBTree(bbs, tmp, level+1, (int)new_elems_right.size(),_epsilon);
 137
 138   }
 139
 140
 141
 142   ~BBTree()
 143   {
 144     if (_left!=0)  delete _left;
 145     if (_right!=0) delete _right;
 146
 147   }
 148
 149
 150   /*! returns in \a elems the list of elements potentially intersecting the bounding box pointed to by \a bb
 151
 152     \param bb pointer to query bounding box
 153     \param elems list of elements (given in 0-indexing that is to say in \b C \b mode) intersecting the bounding box
 154   */
 155
 156   void getIntersectingElems(const double* bb, std::vector<ConnType>& elems) const
 157   {
 158     //  terminal node : return list of elements intersecting bb
 159     if (_terminal)
 160       {
 161         for (int i=0; i<_nbelems; i++)
 162           {
 163             const double* const  bb_ptr=_bb+_elems[i]*2*dim;
 164             bool intersects = true;
 165             for (int idim=0; idim<dim; idim++)
 166               {
 167                 if (bb_ptr[idim*2]-bb[idim*2+1]>-_epsilon|| bb_ptr[idim*2+1]-bb[idim*2]<_epsilon)
 168                   intersects=false;
 169               }
 170             if (intersects)
 171               {
 172                 elems.push_back(_elems[i]);
 173               }
 174           }
 175         return;
 176       }
 177
 178     //non terminal node
 179     double min = bb[(_level%dim)*2];
 180     double max = bb[(_level%dim)*2+1];
 181     if (max < _min_right)
 182       {
 183         _left->getIntersectingElems(bb, elems);
 184         return;
 185       }
 186     if (min> _max_left)
 187       {
 188         _right->getIntersectingElems(bb,elems);
 189         return;
 190       }
 191     _left->getIntersectingElems(bb,elems);
 192     _right->getIntersectingElems(bb,elems);
 193   }
 194
 195
 196   /*! returns in \a elems the list of elements potentially containing the point pointed to by \a xx
 197     \param xx pointer to query point coords
 198     \param elems list of elements (given in 0-indexing) intersecting the bounding box
 199   */
 200
 201   void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const
 202   {
 203     //  terminal node : return list of elements intersecting bb
 204     if (_terminal)
 205       {
 206         for (ConnType 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]-xx[idim]>_epsilon|| bb_ptr[idim*2+1]-xx[idim]<-_epsilon)
 213                   intersects=false;
 214               }
 215             if (intersects)
 216               {
 217                 elems.push_back(_elems[i]);
 218               }
 219           }
 220         return;
 221       }
 222
 223     //non terminal node
 224     if (xx[_level%dim] < _min_right)
 225       {
 226         _left->getElementsAroundPoint(xx, elems);
 227         return;
 228       }
 229     if (xx[_level%dim]> _max_left)
 230       {
 231         _right->getElementsAroundPoint(xx,elems);
 232         return;
 233       }
 234     _left->getElementsAroundPoint(xx,elems);
 235     _right->getElementsAroundPoint(xx,elems);
 236   }
 237
 238
 239
 240   int size()
 241   {
 242     if (_terminal) return _nbelems;
 243     return _left->size()+_right->size();
 244   }
 245 };
 246 #endif