1 // Copyright (C) 2007-2015 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 __BBTREE_TXX__
20 #define __BBTREE_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 bbs pointer to the [xmin1 xmax1 ymin1 ymax1 xmin2 xmax2 ...] array containing the bounding boxes that are to be indexed.
52 \param elems array to the indices of the elements contained in the BBTree
53 \param level level in the BBTree recursive structure
54 \param nbelems nb of elements in the BBTree
55 \param epsilon precision to which points are decided to be coincident. Epsilon can be positive or negative.
56 If \a epsilon is positive the request method will enlarge the computed bounding box (more matching elems return).
57 If negative the given bounding box will be tighten (less matching elems return).
59 Parameters \a elems and \a level are used only by BBTree itself for creating trees recursively. A typical use is therefore :
62 double* bbs= new double[2*2*nbelems];
65 BBTree<2> tree = new BBTree<2>(elems,0,0,nbelems,1e-12);
69 BBTree(const double* bbs, ConnType* elems, int level, ConnType nbelems, double epsilon=1e-12):
70 _left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems),_epsilon(epsilon)
72 if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL)
77 double* nodes=new double [nbelems];
78 _elems.resize(nbelems);
79 for (ConnType i=0; i<nbelems; i++)
88 nodes[i]=bbs[elem*dim*2+(level%dim)*2];
90 if (_terminal) { delete[] nodes; return;}
92 std::nth_element<double*>(nodes, nodes+nbelems/2, nodes+nbelems);
93 double median = *(nodes+nbelems/2);
95 // std:: cout << *median <<std::endl;
97 std::vector<ConnType> new_elems_left;
98 std::vector<ConnType> new_elems_right;
100 new_elems_left.reserve(nbelems/2+1);
101 new_elems_right.reserve(nbelems/2+1);
102 double max_left = -std::numeric_limits<double>::max();
103 double min_right= std::numeric_limits<double>::max();
104 for (int i=0; i<nbelems;i++)
111 double max=bbs[elem*dim*2+(level%dim)*2+1];
112 double min = bbs[elem*dim*2+(level%dim)*2];
116 new_elems_right.push_back(elem);
117 if (min<min_right) min_right = min;
122 new_elems_left.push_back(elem);
123 if (max>max_left) max_left = max;
128 _max_left=max_left+std::abs(_epsilon);
129 _min_right=min_right-std::abs(_epsilon);
132 if(!new_elems_left.empty())
133 tmp=&(new_elems_left[0]);
134 _left=new BBTree(bbs, tmp, level+1, (int)new_elems_left.size(),_epsilon);
136 if(!new_elems_right.empty())
137 tmp=&(new_elems_right[0]);
138 _right=new BBTree(bbs, tmp, level+1, (int)new_elems_right.size(),_epsilon);
146 if (_left!=0) delete _left;
147 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
158 void getIntersectingElems(const double* bb, std::vector<ConnType>& elems) const
160 // terminal node : return list of elements intersecting bb
163 for (int i=0; i<_nbelems; i++)
165 const double* const bb_ptr=_bb+_elems[i]*2*dim;
166 bool intersects = true;
167 for (int idim=0; idim<dim; idim++)
169 if (bb_ptr[idim*2]-bb[idim*2+1]>-_epsilon|| bb_ptr[idim*2+1]-bb[idim*2]<_epsilon)
174 elems.push_back(_elems[i]);
181 double min = bb[(_level%dim)*2];
182 double max = bb[(_level%dim)*2+1];
183 if (max < _min_right)
185 _left->getIntersectingElems(bb, elems);
190 _right->getIntersectingElems(bb,elems);
193 _left->getIntersectingElems(bb,elems);
194 _right->getIntersectingElems(bb,elems);
198 * This method is very close to getIntersectingElems except that it returns number of elems instead of elems themselves.
200 int getNbOfIntersectingElems(const double* bb)
202 // terminal node : return list of elements intersecting bb
206 for (int i=0; i<_nbelems; i++)
208 const double* const bb_ptr=_bb+_elems[i]*2*dim;
209 bool intersects = true;
210 for (int idim=0; idim<dim; idim++)
212 if (bb_ptr[idim*2]-bb[idim*2+1]>-_epsilon|| bb_ptr[idim*2+1]-bb[idim*2]<_epsilon)
221 double min = bb[(_level%dim)*2];
222 double max = bb[(_level%dim)*2+1];
223 if (max < _min_right)
224 return _left->getNbOfIntersectingElems(bb);
226 return _right->getNbOfIntersectingElems(bb);
227 return _left->getNbOfIntersectingElems(bb)+_right->getNbOfIntersectingElems(bb);
231 /*! returns in \a elems the list of elements potentially containing the point pointed to by \a xx
232 \param xx pointer to query point coords
233 \param elems list of elements (given in 0-indexing) intersecting the bounding box
236 void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const
238 // terminal node : return list of elements intersecting bb
241 for (ConnType i=0; i<_nbelems; i++)
243 const double* const bb_ptr=_bb+_elems[i]*2*dim;
244 bool intersects = true;
245 for (int idim=0; idim<dim; idim++)
247 if (bb_ptr[idim*2]-xx[idim]>_epsilon|| bb_ptr[idim*2+1]-xx[idim]<-_epsilon)
252 elems.push_back(_elems[i]);
259 if (xx[_level%dim] < _min_right)
261 _left->getElementsAroundPoint(xx, elems);
264 if (xx[_level%dim]> _max_left)
266 _right->getElementsAroundPoint(xx,elems);
269 _left->getElementsAroundPoint(xx,elems);
270 _right->getElementsAroundPoint(xx,elems);
277 if (_terminal) return _nbelems;
278 return _left->size()+_right->size();