1 // Copyright (C) 2007-2022 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
30 constexpr double BBTREE_DFT_EPSILON = 1e-12;
32 template <int dim, class ConnType = int>
42 typename std::vector<ConnType> _elems;
47 static const int MIN_NB_ELEMS=15;
48 static const int MAX_LEVEL=20;
52 Constructor of the bounding box tree
53 \param bbs pointer to the [xmin1 xmax1 ymin1 ymax1 xmin2 xmax2 ...] array containing the bounding boxes that are to be indexed.
54 \param elems array to the indices of the elements contained in the BBTree
55 \param level level in the BBTree recursive structure
56 \param nbelems nb of elements in the BBTree
57 \param epsilon precision to which points are decided to be coincident. Epsilon can be positive or negative.
58 If \a epsilon is positive the request method will enlarge the computed bounding box (more matching elems return).
59 If negative the given bounding box will be tighten (less matching elems return).
61 Parameters \a elems and \a level are used only by BBTree itself for creating trees recursively. A typical use is therefore :
64 double* bbs= new double[2*2*nbelems];
67 BBTree<2> tree = new BBTree<2>(elems,0,0,nbelems,1e-12);
70 BBTree(const double* bbs, ConnType* elems, int level, ConnType nbelems, double epsilon=BBTREE_DFT_EPSILON):
71 _left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems),_epsilon(epsilon)
73 if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL)
78 double median = std::numeric_limits<double>::max();
80 std::unique_ptr<double[]> nodes( new double [nbelems] );
81 _elems.resize(nbelems);
82 for (ConnType i=0; i<nbelems; i++)
91 nodes[i]=bbs[elem*dim*2+(level%dim)*2];
93 if (_terminal) { return; }
95 std::nth_element<double*>(nodes.get(), nodes.get()+nbelems/2, nodes.get()+nbelems);
96 median = *(nodes.get()+nbelems/2);
98 // std:: cout << *median <<std::endl;
100 std::vector<ConnType> new_elems_left;
101 std::vector<ConnType> new_elems_right;
103 new_elems_left.reserve(nbelems/2+1);
104 new_elems_right.reserve(nbelems/2+1);
105 double max_left = -std::numeric_limits<double>::max();
106 double min_right= std::numeric_limits<double>::max();
107 for (ConnType i=0; i<nbelems;i++)
114 double max=bbs[elem*dim*2+(level%dim)*2+1];
115 double min = bbs[elem*dim*2+(level%dim)*2];
119 new_elems_right.push_back(elem);
120 if (min<min_right) min_right = min;
125 new_elems_left.push_back(elem);
126 if (max>max_left) max_left = max;
131 _max_left=max_left+std::abs(_epsilon);
132 _min_right=min_right-std::abs(_epsilon);
135 if(!new_elems_left.empty())
136 tmp=&(new_elems_left[0]);
137 _left=new BBTree(bbs, tmp, level+1, (ConnType)new_elems_left.size(),_epsilon);
139 if(!new_elems_right.empty())
140 tmp=&(new_elems_right[0]);
141 _right=new BBTree(bbs, tmp, level+1, (ConnType)new_elems_right.size(),_epsilon);
147 if (_left!=0) delete _left;
148 if (_right!=0) delete _right;
152 /*! returns in \a elems the list of elements potentially intersecting the bounding box pointed to by \a bb
154 \param bb pointer to query bounding box
155 \param elems list of elements (given in 0-indexing that is to say in \b C \b mode) intersecting the bounding box
157 void getIntersectingElems(const double* bb, std::vector<ConnType>& elems) const
159 // terminal node : return list of elements intersecting bb
162 for (ConnType i=0; i<_nbelems; i++)
164 const double* const bb_ptr=_bb+_elems[i]*2*dim;
165 bool intersects = true;
166 for (int idim=0; idim<dim; idim++)
168 if (bb_ptr[idim*2]-bb[idim*2+1]>-_epsilon|| bb_ptr[idim*2+1]-bb[idim*2]<_epsilon)
173 elems.push_back(_elems[i]);
180 double min = bb[(_level%dim)*2];
181 double max = bb[(_level%dim)*2+1];
182 if (max < _min_right)
184 _left->getIntersectingElems(bb, elems);
189 _right->getIntersectingElems(bb,elems);
192 _left->getIntersectingElems(bb,elems);
193 _right->getIntersectingElems(bb,elems);
197 * This method is very close to getIntersectingElems except that it returns number of elems instead of elems themselves.
199 ConnType getNbOfIntersectingElems(const double* bb)
201 // terminal node : return list of elements intersecting bb
205 for (ConnType i=0; i<_nbelems; i++)
207 const double* const bb_ptr=_bb+_elems[i]*2*dim;
208 bool intersects = true;
209 for (int idim=0; idim<dim; idim++)
211 if (bb_ptr[idim*2]-bb[idim*2+1]>-_epsilon|| bb_ptr[idim*2+1]-bb[idim*2]<_epsilon)
220 double min = bb[(_level%dim)*2];
221 double max = bb[(_level%dim)*2+1];
222 if (max < _min_right)
223 return _left->getNbOfIntersectingElems(bb);
225 return _right->getNbOfIntersectingElems(bb);
226 return _left->getNbOfIntersectingElems(bb)+_right->getNbOfIntersectingElems(bb);
230 /*! returns in \a elems the list of elements potentially containing the point pointed to by \a xx
231 \param xx pointer to query point coords
232 \param elems list of elements (given in 0-indexing) intersecting the bounding box
234 void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const
236 // terminal node : return list of elements intersecting bb
239 for (ConnType i=0; i<_nbelems; i++)
241 const double* const bb_ptr=_bb+_elems[i]*2*dim;
242 bool intersects = true;
243 for (int idim=0; idim<dim; idim++)
245 if (bb_ptr[idim*2]-xx[idim]>_epsilon|| bb_ptr[idim*2+1]-xx[idim]<-_epsilon)
250 elems.push_back(_elems[i]);
257 if (xx[_level%dim] < _min_right)
259 _left->getElementsAroundPoint(xx, elems);
262 if (xx[_level%dim]> _max_left)
264 _right->getElementsAroundPoint(xx,elems);
267 _left->getElementsAroundPoint(xx,elems);
268 _right->getElementsAroundPoint(xx,elems);
273 if (_terminal) return _nbelems;
274 return _left->size()+_right->size();