1 // Copyright (C) 2007-2013 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.
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(int 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, (int)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, (int)new_elems_right.size(),_epsilon);
133 /*! returns in \a elems the list of elements 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
138 \param elems list of elements (given in 0-indexing) 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);
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 \param xx pointer to query point coords
178 \param elems list of elements (given in 0-indexing) intersecting the bounding box
179 \sa BBTreePts::getElementsAroundPoint2
181 void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const
183 // terminal node : return list of elements intersecting bb
186 for(ConnType i=0;i<_nbelems;i++)
188 const double* const bb_ptr=_pts+_elems[i]*dim;
189 bool intersects = true;
190 for(int idim=0;idim<dim;idim++)
191 if(std::abs(bb_ptr[idim]-xx[idim])>_epsilon)
194 elems.push_back(_elems[i]);
199 if(xx[_level%dim]<_min_right)
201 _left->getElementsAroundPoint(xx,elems);
204 if(xx[_level%dim]>_max_left)
206 _right->getElementsAroundPoint(xx,elems);
209 _left->getElementsAroundPoint(xx,elems);
210 _right->getElementsAroundPoint(xx,elems);
217 return _left->size()+_right->size();