1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #ifndef __BBTREEDST_TXX__
22 #define __BBTREEDST_TXX__
41 std::vector<int> _elems;
45 static const int MIN_NB_ELEMS=15;
46 static const int MAX_LEVEL=20;
48 BBTreeDst(const double* bbs, int* elems, int level, int nbelems):
49 _left(0),_right(0),_level(level),_bb(bbs),_terminal(0),_nbelems(nbelems)
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;
58 fillBBoxTerminal(bbs);
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);
67 std::vector<int> new_elems_left;
68 std::vector<int> new_elems_right;
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++)
81 double max=bbs[elem*dim*2+(level%dim)*2+1];
82 double min = bbs[elem*dim*2+(level%dim)*2];
86 new_elems_right.push_back(elem);
87 if (min<min_right) min_right = min;
91 new_elems_left.push_back(elem);
92 if (max>max_left) max_left = max;
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());
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());
115 void getElemsWhoseMinDistanceToPtSmallerThan(const double *pt, double minOfMaxDstsSq, std::vector<int>& elems) const
119 for(int i=0; i<_nbelems; i++)
121 if(GetMinDistanceFromBBoxToPt(_bb+_elems[i]*2*dim,pt)<minOfMaxDstsSq)
122 elems.push_back(_elems[i]);
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);
137 void getMinDistanceOfMax(const double *pt, double& minOfMaxDstsSq) const
141 if(GetMinDistanceFromBBoxToPt(_terminal,pt)>minOfMaxDstsSq)//min it is not a bug
143 for(int i=0; i<_nbelems; i++)
145 minOfMaxDstsSq=std::min(minOfMaxDstsSq,GetMaxDistanceFromBBoxToPt(_bb+_elems[i]*2*dim,pt));
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);
160 void fillBBoxTerminal(const double* bbs)
162 for(int j=0;j<dim;j++)
164 _terminal[2*j]=std::numeric_limits<double>::max();
165 _terminal[2*j+1]=-std::numeric_limits<double>::max();
167 for(int i=0;i<_nbelems;i++)
169 for(int j=0;j<dim;j++)
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]);
177 static double GetMaxDistanceFromBBoxToPt(const double *bbox, const double *pt)
182 for (int idim=0; idim<dim; idim++)
184 double val1=pt[idim]-bbox[idim*2],val2=pt[idim]-bbox[idim*2+1];
185 double x=std::max(fabs(val1),fabs(val2));
190 else//min>max -> no cells in this
191 return std::numeric_limits<double>::max();
195 static double GetMinDistanceFromBBoxToPt(const double *bbox, const double *pt)
200 for (int idim=0; idim<dim; idim++)
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.)
206 double x=pos==1?val2:val1;
212 else//min>max -> no cells in this
213 return std::numeric_limits<double>::max();