1 // Copyright (C) 2007-2023 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 #ifndef __BBTREEPTS_TXX__
20 #define __BBTREEPTS_TXX__
29 template <int dim, class ConnType = int>
40 typename std::vector<ConnType> _elems;
45 static const int MIN_NB_ELEMS=15;
46 static const int MAX_LEVEL=20;
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.
57 Parameters \a elems and \a level are used only by BBTreePts itself for creating trees recursively. A typical use is therefore :
60 double* pts= new double[dim*nbelems];
63 BBTreePts<2> tree = new BBTreePts<2>(elems,0,0,nbelems,1e-12);
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))
69 double *nodes=new double[nbelems];
70 _elems.resize(nbelems);
71 for (ConnType i=0;i<nbelems;i++)
80 nodes[i]=pts[elem*dim+(level%dim)];
82 if(_terminal) { delete[] nodes; return; }
84 std::nth_element<double*>(nodes, nodes+nbelems/2, nodes+nbelems);
85 double median=*(nodes+nbelems/2);
87 std::vector<ConnType> new_elems_left,new_elems_right;
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++)
100 double mx=pts[elem*dim+(level%dim)];
103 new_elems_right.push_back(elem);
104 if(mx<min_right) min_right=mx;
108 new_elems_left.push_back(elem);
109 if(mx>max_left) max_left=mx;
112 _max_left=max_left+_epsilon;
113 _min_right=min_right-_epsilon;
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);
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);
133 /*! returns in \a elem an element potentially containing the point pointed to by \a xx
134 Contrary to BBTreePts::getElementsAroundPoint the norm 2 is used here.
136 \param[in] xx pointer to query point coords
137 \param[in] threshold detection precision
138 \param[out] elem the element intersecting the bounding box
139 \sa BBTreePts::getElementsAroundPoint
141 double getElementsAroundPoint2(const double *xx, double threshold, ConnType& elem) const
143 // terminal node : return list of elements intersecting bb
146 double ret=std::numeric_limits<double>::max();
147 for(ConnType i=0;i<_nbelems;i++)
149 const double* const bb_ptr=_pts+_elems[i]*dim;
151 for(int idim=0;idim<dim;idim++)
152 tmp+=(bb_ptr[idim]-xx[idim])*(bb_ptr[idim]-xx[idim]);
156 { ret=tmp; elem=_elems[i]; }
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);
171 { elem=eleml; return retl; }
173 { elem=elemr; return retr; }
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 ! ***
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
183 void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const
185 // terminal node : return list of elements intersecting bb
188 for(ConnType i=0;i<_nbelems;i++)
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);
195 elems.push_back(_elems[i]);
200 if(xx[_level%dim]<_min_right)
202 _left->getElementsAroundPoint(xx,elems);
205 if(xx[_level%dim]>_max_left)
207 _right->getElementsAroundPoint(xx,elems);
210 _left->getElementsAroundPoint(xx,elems);
211 _right->getElementsAroundPoint(xx,elems);
214 ConnType size() const
218 return _left->size()+_right->size();