src/INTERP_KERNEL/BBTree.txx

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