src/INTERP_KERNEL/BBTreeDst.txx

   1 // Copyright (C) 2007-2016  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 // Author : Anthony Geay (CEA/DEN)
  20
  21 #ifndef __BBTREEDST_TXX__
  22 #define __BBTREEDST_TXX__
  23
  24 #include <vector>
  25 #include <algorithm>
  26
  27 #include <iostream>
  28 #include <limits>
  29 #include <cmath>
  30
  31 template <int dim>
  32 class BBTreeDst
  33 {
  34 private:
  35   BBTreeDst* _left;
  36   BBTreeDst* _right;
  37   int _level;
  38   double _max_left;
  39   double _min_right;
  40   const double *_bb;
  41   std::vector<int> _elems;
  42   double  *_terminal;
  43   int _nbelems;
  44
  45   static const int MIN_NB_ELEMS=15;
  46   static const int MAX_LEVEL=20;
  47 public:
  48   BBTreeDst(const double* bbs, int* elems, int level, int nbelems):
  49     _left(0),_right(0),_level(level),_bb(bbs),_terminal(0),_nbelems(nbelems)
  50   {
  51     if((nbelems < MIN_NB_ELEMS || level> MAX_LEVEL))
  52       _terminal=new double[2*dim];
  53     _elems.resize(nbelems);
  54     for (int i=0; i<nbelems; i++)
  55       _elems[i]=elems?elems[i]:i;
  56     if(_terminal)
  57       {
  58         fillBBoxTerminal(bbs);
  59         return ;
  60       }
  61     double *nodes=new double[nbelems];
  62     for (int i=0; i<nbelems; i++)
  63       nodes[i]=bbs[_elems[i]*dim*2+(level%dim)*2];
  64     std::nth_element<double*>(nodes, nodes+nbelems/2, nodes+nbelems);
  65     double median = *(nodes+nbelems/2);
  66     delete [] nodes;
  67     std::vector<int> new_elems_left;
  68     std::vector<int> new_elems_right;
  69
  70     new_elems_left.reserve(nbelems/2+1);
  71     new_elems_right.reserve(nbelems/2+1);
  72     double max_left = -std::numeric_limits<double>::max();
  73     double min_right=  std::numeric_limits<double>::max();
  74     for(int i=0; i<nbelems;i++)
  75       {
  76         int elem;
  77         if (elems!=0)
  78           elem= elems[i];
  79         else
  80           elem=i;
  81         double max=bbs[elem*dim*2+(level%dim)*2+1];
  82         double min = bbs[elem*dim*2+(level%dim)*2];
  83
  84         if (min>median)
  85           {
  86             new_elems_right.push_back(elem);
  87             if (min<min_right) min_right = min;
  88           }
  89         else
  90           {
  91             new_elems_left.push_back(elem);
  92             if (max>max_left) max_left = max;
  93           }
  94       }
  95     _max_left=max_left;
  96     _min_right=min_right;
  97     int *tmp;
  98     tmp=0;
  99     if(!new_elems_left.empty())
 100       tmp=&(new_elems_left[0]);
 101     _left=new BBTreeDst(bbs, tmp, level+1, (int)new_elems_left.size());
 102     tmp=0;
 103     if(!new_elems_right.empty())
 104       tmp=&(new_elems_right[0]);
 105     _right=new BBTreeDst(bbs, tmp, level+1, (int)new_elems_right.size());
 106   }
 107
 108   ~BBTreeDst()
 109   {
 110     delete _left;
 111     delete _right;
 112     delete [] _terminal;
 113   }
 114
 115   void getElemsWhoseMinDistanceToPtSmallerThan(const double *pt, double minOfMaxDstsSq, std::vector<int>& elems) const
 116   {
 117     if(_terminal)
 118       {
 119         for(int i=0; i<_nbelems; i++)
 120           {
 121             if(GetMinDistanceFromBBoxToPt(_bb+_elems[i]*2*dim,pt)<minOfMaxDstsSq)
 122               elems.push_back(_elems[i]);
 123           }
 124       }
 125     else
 126       {
 127         double minOfMaxDsts=sqrt(minOfMaxDstsSq);
 128         if(_min_right-pt[_level%dim]>minOfMaxDsts)
 129           { _left->getElemsWhoseMinDistanceToPtSmallerThan(pt,minOfMaxDstsSq,elems); return ; }
 130         if(pt[_level%dim]-_max_left>minOfMaxDsts)
 131           { _right->getElemsWhoseMinDistanceToPtSmallerThan(pt,minOfMaxDstsSq,elems); return ; }
 132         _left->getElemsWhoseMinDistanceToPtSmallerThan(pt,minOfMaxDstsSq,elems);
 133         _right->getElemsWhoseMinDistanceToPtSmallerThan(pt,minOfMaxDstsSq,elems);
 134       }
 135   }
 136
 137   /** Get the minimal (square) distance between a point and all the available bounding boxes in the tree.
 138     The (square) distance to a bbox is the true geometric distance between the point and a face
 139     (or an edge, or a corner) of the bbox.
 140   */
 141   void getMinDistanceOfMax(const double *pt, double& minOfMaxDstsSq) const
 142   {
 143     if(_terminal)
 144       {
 145         if(GetMinDistanceFromBBoxToPt(_terminal,pt)>minOfMaxDstsSq)//min it is not a bug
 146           return ;
 147         for(int i=0; i<_nbelems; i++)
 148           {
 149             minOfMaxDstsSq=std::min(minOfMaxDstsSq,GetMaxDistanceFromBBoxToPt(_bb+_elems[i]*2*dim,pt));
 150           }
 151       }
 152     else
 153       {
 154         double minOfMaxDsts=sqrt(minOfMaxDstsSq);
 155         if(_min_right-pt[_level%dim]>minOfMaxDsts)
 156           { _left->getMinDistanceOfMax(pt,minOfMaxDstsSq); return ; }
 157         if(pt[_level%dim]-_max_left>minOfMaxDsts)
 158           { _right->getMinDistanceOfMax(pt,minOfMaxDstsSq); return ; }
 159         _left->getMinDistanceOfMax(pt,minOfMaxDstsSq);
 160         _right->getMinDistanceOfMax(pt,minOfMaxDstsSq);
 161       }
 162   }
 163
 164   void fillBBoxTerminal(const double* bbs)
 165   {
 166     for(int j=0;j<dim;j++)
 167       {
 168         _terminal[2*j]=std::numeric_limits<double>::max();
 169         _terminal[2*j+1]=-std::numeric_limits<double>::max();
 170       }
 171     for(int i=0;i<_nbelems;i++)
 172       {
 173         for(int j=0;j<dim;j++)
 174           {
 175             _terminal[2*j]=std::min(_terminal[2*j],bbs[2*dim*_elems[i]+2*j]);
 176             _terminal[2*j+1]=std::max(_terminal[2*j+1],bbs[2*dim*_elems[i]+2*j+1]);
 177           }
 178       }
 179   }
 180
 181   static double GetMaxDistanceFromBBoxToPt(const double *bbox, const double *pt)
 182   {
 183     if(bbox[0]<=bbox[1])
 184       {
 185         double zeRes=0.;
 186         for (int idim=0; idim<dim; idim++)
 187           {
 188             double val1=pt[idim]-bbox[idim*2],val2=pt[idim]-bbox[idim*2+1];
 189             double x=std::max(fabs(val1),fabs(val2));
 190             zeRes+=x*x;
 191           }
 192         return zeRes;
 193       }
 194     else//min>max -> no cells in this
 195       return std::numeric_limits<double>::max();
 196
 197   }
 198
 199   static double GetMinDistanceFromBBoxToPt(const double *bbox, const double *pt)
 200   {
 201     if(bbox[0]<=bbox[1])
 202       {
 203         double zeRes=0.;
 204         for (int idim=0; idim<dim; idim++)
 205           {
 206             double val1=pt[idim]-bbox[idim*2],val2=pt[idim]-bbox[idim*2+1];
 207             char pos=(( (0.<val1)-(val1<0.) )+( (0.<val2)-(val2<0.) ))/2;// sign(val) = (0.<val)-(val<0.)
 208             if(pos!=0)
 209               {
 210                 double x=pos==1?val2:val1;
 211                 zeRes+=x*x;
 212               }
 213           }
 214         return zeRes;
 215       }
 216     else//min>max -> no cells in this
 217       return std::numeric_limits<double>::max();
 218   }
 219 };
 220
 221 #endif