src/INTERP_KERNEL/BBTreePts.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 __BBTREEPTS_TXX__
  20 #define __BBTREEPTS_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 BBTreePts
  31 {
  32
  33 private:
  34   BBTreePts* _left;
  35   BBTreePts* _right;
  36   int _level;
  37   double _max_left;
  38   double _min_right;
  39   const double *_pts;
  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 [in] pts pointer to the array containing the points that are to be indexed.
  52     \param [in] elems array to the indices of the elements contained in the BBTreePts
  53     \param [in] level level in the BBTreePts recursive structure
  54     \param [in] nbelems nb of elements in the BBTreePts
  55     \param [in] epsilon precision to which points are decided to be coincident. Contrary to BBTree, the absolute epsilon is computed. So the internal epsilon is always positive.
  56
  57     Parameters \a elems and \a level are used only by BBTreePts itself for creating trees recursively. A typical use is therefore :
  58     \code
  59     int nbelems=...
  60     double* pts= new double[dim*nbelems];
  61     // filling pts ...
  62     ...
  63     BBTreePts<2> tree = new BBTreePts<2>(elems,0,0,nbelems,1e-12);
  64     \endcode
  65   */
  66   BBTreePts(const double *pts, const ConnType *elems, int level, ConnType nbelems, double epsilon=1e-12):
  67     _left(0),_right(0),_level(level),_pts(pts),_terminal(nbelems < MIN_NB_ELEMS || level> MAX_LEVEL),_nbelems(nbelems),_epsilon(std::abs(epsilon))
  68   {
  69     double *nodes=new double[nbelems];
  70     _elems.resize(nbelems);
  71     for (ConnType i=0;i<nbelems;i++)
  72       {
  73         ConnType elem;
  74         if (elems!=0)
  75           elem= elems[i];
  76         else
  77           elem=i;
  78
  79         _elems[i]=elem;
  80         nodes[i]=pts[elem*dim+(level%dim)];
  81       }
  82     if(_terminal) { delete[] nodes; return; }
  83     //
  84     std::nth_element<double*>(nodes, nodes+nbelems/2, nodes+nbelems);
  85     double median=*(nodes+nbelems/2);
  86     delete [] nodes;
  87     std::vector<ConnType> new_elems_left,new_elems_right;
  88
  89     new_elems_left.reserve(nbelems/2+1);
  90     new_elems_right.reserve(nbelems/2+1);
  91     double max_left = -std::numeric_limits<double>::max();
  92     double min_right=  std::numeric_limits<double>::max();
  93     for(int i=0;i<nbelems;i++)
  94       {
  95         int elem;
  96         if(elems!=0)
  97           elem= elems[i];
  98         else
  99           elem=i;
 100         double mx=pts[elem*dim+(level%dim)];
 101         if(mx>median)
 102           {
 103             new_elems_right.push_back(elem);
 104             if(mx<min_right) min_right=mx;
 105           }
 106         else
 107           {
 108             new_elems_left.push_back(elem);
 109             if(mx>max_left) max_left=mx;
 110           }
 111       }
 112     _max_left=max_left+_epsilon;
 113     _min_right=min_right-_epsilon;
 114     ConnType *tmp;
 115     tmp=0;
 116     if(!new_elems_left.empty())
 117       tmp=&(new_elems_left[0]);
 118     _left=new BBTreePts(pts, tmp, level+1, (int)new_elems_left.size(),_epsilon);
 119     tmp=0;
 120     if(!new_elems_right.empty())
 121       tmp=&(new_elems_right[0]);
 122     _right=new BBTreePts(pts, tmp, level+1, (int)new_elems_right.size(),_epsilon);
 123   }
 124
 125
 126
 127   ~BBTreePts()
 128   {
 129     delete _left;
 130     delete _right;
 131   }
 132
 133   /*! returns in \a elems the list of elements potentially containing the point pointed to by \a xx
 134       Contrary to BBTreePts::getElementsAroundPoint the norm 2 is used here.
 135
 136     \param [in] xx pointer to query point coords
 137     \param [in] threshold
 138     \param elems list of elements (given in 0-indexing) intersecting the bounding box
 139     \sa BBTreePts::getElementsAroundPoint
 140   */
 141   double getElementsAroundPoint2(const double *xx, double threshold, ConnType& elem) const
 142   {
 143     //  terminal node : return list of elements intersecting bb
 144     if(_terminal)
 145       {
 146         double ret=std::numeric_limits<double>::max();
 147         for(ConnType i=0;i<_nbelems;i++)
 148           {
 149             const double* const bb_ptr=_pts+_elems[i]*dim;
 150             double tmp=0.;
 151             for(int idim=0;idim<dim;idim++)
 152               tmp+=(bb_ptr[idim]-xx[idim])*(bb_ptr[idim]-xx[idim]);
 153             if(tmp<threshold)
 154               {
 155                 if(tmp<ret)
 156                   { ret=tmp; elem=_elems[i]; }
 157               }
 158           }
 159         return ret;
 160       }
 161     //non terminal node
 162     double s=sqrt(threshold*dim);
 163     if(xx[_level%dim]+s<_min_right)
 164       return _left->getElementsAroundPoint2(xx,threshold,elem);
 165     if(xx[_level%dim]-s>_max_left)
 166       return _right->getElementsAroundPoint2(xx,threshold,elem);
 167     int eleml,elemr;
 168     double retl=_left->getElementsAroundPoint2(xx,threshold,eleml);
 169     double retr=_right->getElementsAroundPoint2(xx,threshold,elemr);
 170     if(retl<retr)
 171       { elem=eleml; return retl; }
 172     else
 173       { elem=elemr; return retr; }
 174   }
 175
 176   /*! returns in \a elems the list of elements potentially containing the point pointed to by \a xx
 177     \param xx pointer to query point coords
 178     \param elems list of elements (given in 0-indexing) intersecting the bounding box
 179     \sa BBTreePts::getElementsAroundPoint2
 180   */
 181   void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const
 182   {
 183     //  terminal node : return list of elements intersecting bb
 184     if(_terminal)
 185       {
 186         for(ConnType i=0;i<_nbelems;i++)
 187           {
 188             const double* const bb_ptr=_pts+_elems[i]*dim;
 189             bool intersects = true;
 190             for(int idim=0;idim<dim;idim++)
 191               if(std::abs(bb_ptr[idim]-xx[idim])>_epsilon)
 192                 intersects=false;
 193             if(intersects)
 194               elems.push_back(_elems[i]);
 195           }
 196         return;
 197       }
 198     //non terminal node
 199     if(xx[_level%dim]<_min_right)
 200       {
 201         _left->getElementsAroundPoint(xx,elems);
 202         return;
 203       }
 204     if(xx[_level%dim]>_max_left)
 205       {
 206         _right->getElementsAroundPoint(xx,elems);
 207         return;
 208       }
 209     _left->getElementsAroundPoint(xx,elems);
 210     _right->getElementsAroundPoint(xx,elems);
 211   }
 212
 213   int size() const
 214   {
 215     if(_terminal)
 216       return _nbelems;
 217     return _left->size()+_right->size();
 218   }
 219 };
 220
 221 #endif