src/INTERP_KERNEL/BBTreeDst.txx

   1 // Copyright (C) 2007-2015  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   void getMinDistanceOfMax(const double *pt, double& minOfMaxDstsSq) const
 138   {
 139     if(_terminal)
 140       {
 141         if(GetMinDistanceFromBBoxToPt(_terminal,pt)>minOfMaxDstsSq)//min it is not a bug
 142           return ;
 143         for(int i=0; i<_nbelems; i++)
 144           {
 145             minOfMaxDstsSq=std::min(minOfMaxDstsSq,GetMaxDistanceFromBBoxToPt(_bb+_elems[i]*2*dim,pt));
 146           }
 147       }
 148     else
 149       {
 150         double minOfMaxDsts=sqrt(minOfMaxDstsSq);
 151         if(_min_right-pt[_level%dim]>minOfMaxDsts)
 152           { _left->getMinDistanceOfMax(pt,minOfMaxDstsSq); return ; }
 153         if(pt[_level%dim]-_max_left>minOfMaxDsts)
 154           { _right->getMinDistanceOfMax(pt,minOfMaxDstsSq); return ; }
 155         _left->getMinDistanceOfMax(pt,minOfMaxDstsSq);
 156         _right->getMinDistanceOfMax(pt,minOfMaxDstsSq);
 157       }
 158   }
 159
 160   void fillBBoxTerminal(const double* bbs)
 161   {
 162     for(int j=0;j<dim;j++)
 163       {
 164         _terminal[2*j]=std::numeric_limits<double>::max();
 165         _terminal[2*j+1]=-std::numeric_limits<double>::max();
 166       }
 167     for(int i=0;i<_nbelems;i++)
 168       {
 169         for(int j=0;j<dim;j++)
 170           {
 171             _terminal[2*j]=std::min(_terminal[2*j],bbs[2*dim*_elems[i]+2*j]);
 172             _terminal[2*j+1]=std::max(_terminal[2*j+1],bbs[2*dim*_elems[i]+2*j+1]);
 173           }
 174       }
 175   }
 176
 177   static double GetMaxDistanceFromBBoxToPt(const double *bbox, const double *pt)
 178   {
 179     if(bbox[0]<=bbox[1])
 180       {
 181         double zeRes=0.;
 182         for (int idim=0; idim<dim; idim++)
 183           {
 184             double val1=pt[idim]-bbox[idim*2],val2=pt[idim]-bbox[idim*2+1];
 185             double x=std::max(fabs(val1),fabs(val2));
 186             zeRes+=x*x;
 187           }
 188         return zeRes;
 189       }
 190     else//min>max -> no cells in this
 191       return std::numeric_limits<double>::max();
 192
 193   }
 194
 195   static double GetMinDistanceFromBBoxToPt(const double *bbox, const double *pt)
 196   {
 197     if(bbox[0]<=bbox[1])
 198       {
 199         double zeRes=0.;
 200         for (int idim=0; idim<dim; idim++)
 201           {
 202             double val1=pt[idim]-bbox[idim*2],val2=pt[idim]-bbox[idim*2+1];
 203             char pos=(( (0.<val1)-(val1<0.) )+( (0.<val2)-(val2<0.) ))/2;// sign(val) = (0.<val)-(val<0.)
 204             if(pos!=0)
 205               {
 206                 double x=pos==1?val2:val1;
 207                 zeRes+=x*x;
 208               }
 209           }
 210         return zeRes;
 211       }
 212     else//min>max -> no cells in this
 213       return std::numeric_limits<double>::max();
 214   }
 215 };
 216
 217 #endif