src/INTERP_KERNEL/BBTreePts.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 #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(ConnType i=0;i<nbelems;i++)
  94       {
  95         ConnType 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, (ConnType)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, (ConnType)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     ConnType 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    * ** Infinite norm is used here not norm 2 ! ***
 178    *
 179    *  \param xx pointer to query point coords
 180    *  \param elems list of elements (given in 0-indexing) intersecting the bounding box
 181    * \sa BBTreePts::getElementsAroundPoint2
 182    */
 183   void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const
 184   {
 185     //  terminal node : return list of elements intersecting bb
 186     if(_terminal)
 187       {
 188         for(ConnType i=0;i<_nbelems;i++)
 189           {
 190             const double* const bb_ptr=_pts+_elems[i]*dim;
 191             bool intersects = true;
 192             for(int idim=0;idim<dim;idim++)
 193               intersects=intersects && (std::abs(bb_ptr[idim]-xx[idim])<=_epsilon);
 194             if(intersects)
 195               elems.push_back(_elems[i]);
 196           }
 197         return;
 198       }
 199     //non terminal node
 200     if(xx[_level%dim]<_min_right)
 201       {
 202         _left->getElementsAroundPoint(xx,elems);
 203         return;
 204       }
 205     if(xx[_level%dim]>_max_left)
 206       {
 207         _right->getElementsAroundPoint(xx,elems);
 208         return;
 209       }
 210     _left->getElementsAroundPoint(xx,elems);
 211     _right->getElementsAroundPoint(xx,elems);
 212   }
 213
 214   ConnType size() const
 215   {
 216     if(_terminal)
 217       return _nbelems;
 218     return _left->size()+_right->size();
 219   }
 220 };
 221
 222 #endif