1 // Copyright (C) 2007-2024 CEA, EDF
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<mcIdType> _elems;
45 static const int MIN_NB_ELEMS=15;
46 static const int MAX_LEVEL=20;
48 BBTreeDst(const double* bbs, mcIdType* elems, int level, mcIdType 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 (mcIdType i=0; i<nbelems; i++)
55 _elems[i]=elems?elems[i]:i;
58 fillBBoxTerminal(bbs);
61 double *nodes=new double[nbelems];
62 for (mcIdType 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<mcIdType> new_elems_left;
68 std::vector<mcIdType> 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(mcIdType 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, ToIdType(new_elems_left.size()));
103 if(!new_elems_right.empty())
104 tmp=&(new_elems_right[0]);
105 _right=new BBTreeDst(bbs, tmp, level+1, ToIdType(new_elems_right.size()));
115 void getElemsWhoseMinDistanceToPtSmallerThan(const double *pt, double minOfMaxDstsSq, std::vector<mcIdType>& elems) const
119 for(mcIdType 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 /** 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.
141 void getMinDistanceOfMax(const double *pt, double& minOfMaxDstsSq) const
145 if(GetMinDistanceFromBBoxToPt(_terminal,pt)>minOfMaxDstsSq)//min it is not a bug
147 for(mcIdType i=0; i<_nbelems; i++)
149 minOfMaxDstsSq=std::min(minOfMaxDstsSq,GetMaxDistanceFromBBoxToPt(_bb+_elems[i]*2*dim,pt));
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);
164 void fillBBoxTerminal(const double* bbs)
166 for(int j=0;j<dim;j++)
168 _terminal[2*j]=std::numeric_limits<double>::max();
169 _terminal[2*j+1]=-std::numeric_limits<double>::max();
171 for(mcIdType i=0;i<_nbelems;i++)
173 for(int j=0;j<dim;j++)
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]);
181 static double GetMaxDistanceFromBBoxToPt(const double *bbox, const double *pt)
186 for (int idim=0; idim<dim; idim++)
188 double val1=pt[idim]-bbox[idim*2],val2=pt[idim]-bbox[idim*2+1];
189 double x=std::max(fabs(val1),fabs(val2));
194 else//min>max -> no cells in this
195 return std::numeric_limits<double>::max();
199 static double GetMinDistanceFromBBoxToPt(const double *bbox, const double *pt)
204 for (int idim=0; idim<dim; idim++)
206 double val1=pt[idim]-bbox[idim*2],val2=pt[idim]-bbox[idim*2+1];
207 char pos=static_cast<char>((( (0.<val1)-(val1<0.) )+( (0.<val2)-(val2<0.) ))/2);// sign(val) = (0.<val)-(val<0.)
210 double x=pos==1?val2:val1;
216 else//min>max -> no cells in this
217 return std::numeric_limits<double>::max();