\param elems array to the indices of the elements contained in the BBTree
\param level level in the BBTree recursive structure
\param nbelems nb of elements in the BBTree
- \param epsilon precision to which points are decided to be coincident
+ \param epsilon precision to which points are decided to be coincident. Epsilon can be positive or negative.
+ If \a epsilon is positive the request method will enlarge the computed bounding box (more matching elems return).
+ If negative the given bounding box will be tighten (less matching elems return).
Parameters \a elems and \a level are used only by BBTree itself for creating trees recursively. A typical use is therefore :
\code
\endcode
*/
- BBTree(const double* bbs, ConnType* elems, int level, ConnType nbelems, double epsilon=1E-12):
+ BBTree(const double* bbs, ConnType* elems, int level, ConnType nbelems, double epsilon=1e-12):
_left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems),_epsilon(epsilon)
{
if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL)
--- /dev/null
+// Copyright (C) 2007-2012 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+#ifndef __BBTREEPTS_TXX__
+#define __BBTREEPTS_TXX__
+
+#include <vector>
+#include <algorithm>
+
+#include <iostream>
+#include <limits>
+#include <cmath>
+
+template <int dim, class ConnType = int>
+class BBTreePts
+{
+
+private:
+ BBTreePts* _left;
+ BBTreePts* _right;
+ int _level;
+ double _max_left;
+ double _min_right;
+ const double *_pts;
+ typename std::vector<ConnType> _elems;
+ bool _terminal;
+ ConnType _nbelems;
+ double _epsilon;
+
+ static const int MIN_NB_ELEMS=15;
+ static const int MAX_LEVEL=20;
+public:
+
+ /*!
+ Constructor of the bounding box tree
+ \param [in] pts pointer to the array containing the points that are to be indexed.
+ \param [in] elems array to the indices of the elements contained in the BBTreePts
+ \param [in] level level in the BBTreePts recursive structure
+ \param [in] nbelems nb of elements in the BBTreePts
+ \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.
+
+ Parameters \a elems and \a level are used only by BBTreePts itself for creating trees recursively. A typical use is therefore :
+ \code
+ int nbelems=...
+ double* pts= new double[dim*nbelems];
+ // filling pts ...
+ ...
+ BBTreePts<2> tree = new BBTreePts<2>(elems,0,0,nbelems,1e-12);
+ \endcode
+ */
+ BBTreePts(const double *pts, const ConnType *elems, int level, ConnType nbelems, double epsilon=1e-12):
+ _left(0),_right(0),_level(level),_pts(pts),_terminal(nbelems < MIN_NB_ELEMS || level> MAX_LEVEL),_nbelems(nbelems),_epsilon(std::abs(epsilon))
+ {
+ double *nodes=new double[nbelems];
+ _elems.resize(nbelems);
+ for (ConnType i=0;i<nbelems;i++)
+ {
+ ConnType elem;
+ if (elems!=0)
+ elem= elems[i];
+ else
+ elem=i;
+
+ _elems[i]=elem;
+ nodes[i]=pts[elem*dim+(level%dim)];
+ }
+ if(_terminal) { delete[] nodes; return; }
+ //
+ std::nth_element<double*>(nodes, nodes+nbelems/2, nodes+nbelems);
+ double median=*(nodes+nbelems/2);
+ delete [] nodes;
+ std::vector<ConnType> new_elems_left,new_elems_right;
+
+ new_elems_left.reserve(nbelems/2+1);
+ new_elems_right.reserve(nbelems/2+1);
+ double max_left = -std::numeric_limits<double>::max();
+ double min_right= std::numeric_limits<double>::max();
+ for(int i=0;i<nbelems;i++)
+ {
+ int elem;
+ if(elems!=0)
+ elem= elems[i];
+ else
+ elem=i;
+ double mx=pts[elem*dim+(level%dim)];
+ if(mx>median)
+ {
+ new_elems_right.push_back(elem);
+ if(mx<min_right) min_right=mx;
+ }
+ else
+ {
+ new_elems_left.push_back(elem);
+ if(mx>max_left) max_left=mx;
+ }
+ }
+ _max_left=max_left+_epsilon;
+ _min_right=min_right-_epsilon;
+ ConnType *tmp;
+ tmp=0;
+ if(!new_elems_left.empty())
+ tmp=&(new_elems_left[0]);
+ _left=new BBTreePts(pts, tmp, level+1, (int)new_elems_left.size(),_epsilon);
+ tmp=0;
+ if(!new_elems_right.empty())
+ tmp=&(new_elems_right[0]);
+ _right=new BBTreePts(pts, tmp, level+1, (int)new_elems_right.size(),_epsilon);
+ }
+
+
+
+ ~BBTreePts()
+ {
+ delete _left;
+ delete _right;
+ }
+
+ /*! returns in \a elems the list of elements potentially containing the point pointed to by \a xx
+ Contrary to BBTreePts::getElementsAroundPoint the norm 2 is used here.
+
+ \param [in] xx pointer to query point coords
+ \param [in] threshold
+ \param elems list of elements (given in 0-indexing) intersecting the bounding box
+ \sa BBTreePts::getElementsAroundPoint
+ */
+ double getElementsAroundPoint2(const double *xx, double threshold, ConnType& elem) const
+ {
+ // terminal node : return list of elements intersecting bb
+ if(_terminal)
+ {
+ double ret=std::numeric_limits<double>::max();
+ for(ConnType i=0;i<_nbelems;i++)
+ {
+ const double* const bb_ptr=_pts+_elems[i]*dim;
+ double tmp=0.;
+ for(int idim=0;idim<dim;idim++)
+ tmp+=(bb_ptr[idim]-xx[idim])*(bb_ptr[idim]-xx[idim]);
+ if(tmp<threshold)
+ {
+ if(tmp<ret)
+ { ret=tmp; elem=_elems[i]; }
+ }
+ }
+ return ret;
+ }
+ //non terminal node
+ double s=sqrt(threshold*dim);
+ if(xx[_level%dim]+s<_min_right)
+ return _left->getElementsAroundPoint2(xx,threshold,elem);
+ if(xx[_level%dim]-s>_max_left)
+ return _right->getElementsAroundPoint2(xx,threshold,elem);
+ int eleml,elemr;
+ double retl=_left->getElementsAroundPoint2(xx,threshold,eleml);
+ double retr=_right->getElementsAroundPoint2(xx,threshold,elemr);
+ if(retl<retr)
+ { elem=eleml; return retl; }
+ else
+ { elem=elemr; return retr; }
+ }
+
+ /*! returns in \a elems the list of elements potentially containing the point pointed to by \a xx
+ \param xx pointer to query point coords
+ \param elems list of elements (given in 0-indexing) intersecting the bounding box
+ \sa BBTreePts::getElementsAroundPoint2
+ */
+ void getElementsAroundPoint(const double* xx, std::vector<ConnType>& elems) const
+ {
+ // terminal node : return list of elements intersecting bb
+ if(_terminal)
+ {
+ for(ConnType i=0;i<_nbelems;i++)
+ {
+ const double* const bb_ptr=_pts+_elems[i]*dim;
+ bool intersects = true;
+ for(int idim=0;idim<dim;idim++)
+ if(std::abs(bb_ptr[idim]-xx[idim])>_epsilon)
+ intersects=false;
+ if(intersects)
+ elems.push_back(_elems[i]);
+ }
+ return;
+ }
+ //non terminal node
+ if(xx[_level%dim]<_min_right)
+ {
+ _left->getElementsAroundPoint(xx,elems);
+ return;
+ }
+ if(xx[_level%dim]>_max_left)
+ {
+ _right->getElementsAroundPoint(xx,elems);
+ return;
+ }
+ _left->getElementsAroundPoint(xx,elems);
+ _right->getElementsAroundPoint(xx,elems);
+ }
+
+ int size() const
+ {
+ if(_terminal)
+ return _nbelems;
+ return _left->size()+_right->size();
+ }
+};
+
+#endif
salomeinclude_HEADERS = \
BBTree.txx \
+BBTreePts.txx \
BoundingBox.hxx \
CellModel.hxx \
ConvexIntersector.hxx \
}
template<int SPACEDIM>
-int DataArrayDouble::ComputeBasicClosestTupleIdAlg(const std::vector<int>& elems, const double *thisPt, const double *zePt)
-{
- double val=std::numeric_limits<double>::max();
- int ret=-1;
- for(std::vector<int>::const_iterator it=elems.begin();it!=elems.end();it++)
- {
- double tmp=0.;
- for(int j=0;j<SPACEDIM;j++) tmp+=thisPt[SPACEDIM*(*it)+j]-zePt[j];
- if(tmp<val)
- { val=tmp; ret=*it; }
- }
- return ret;
-}
-
-template<int SPACEDIM>
-void DataArrayDouble::FindClosestTupleIdAlg(const BBTree<SPACEDIM,int>& myTree, double dist, const double *pos, int nbOfTuples, const double *thisPt, int thisNbOfTuples, int *res)
+void DataArrayDouble::FindClosestTupleIdAlg(const BBTreePts<SPACEDIM,int>& myTree, double dist, const double *pos, int nbOfTuples, const double *thisPt, int thisNbOfTuples, int *res)
{
double distOpt(dist);
const double *p(pos);
int *r(res);
- double bbox[2*SPACEDIM];
for(int i=0;i<nbOfTuples;i++,p+=SPACEDIM,r++)
{
- double failVal(distOpt),okVal(std::numeric_limits<double>::max());
int nbOfTurn=15;
- while(nbOfTurn>0)
- {
- for(int j=0;j<SPACEDIM;j++) { bbox[2*j]=p[j]-distOpt; bbox[2*j]=p[j]+distOpt; }
- std::vector<int> elems;
- myTree.getIntersectingElems(bbox,elems); nbOfTurn--;
- if(elems.empty())
- { failVal=distOpt; distOpt=okVal==std::numeric_limits<double>::max()?2*distOpt:(okVal+failVal)/2.; continue; }
- if( elems.size()<=15 || nbOfTurn>0)
- { *r=ComputeBasicClosestTupleIdAlg<SPACEDIM>(elems,thisPt,p); break; }
+ while(true)
+ {
+ int elem=-1;
+ double ret=myTree.getElementsAroundPoint2(p,distOpt,elem); nbOfTurn--;
+ if(ret!=std::numeric_limits<double>::max())
+ {
+ distOpt=ret>1e-4?ret:dist;
+ *r=elem;
+ break;
+ }
else
- { okVal=distOpt; distOpt=(distOpt+failVal)/2.; continue; }
+ { distOpt=2*distOpt; continue; }
}
}
}
return ret.retn();
}
+/*!
+ * This methods returns the minimal distance between the two set of points \a this and \a other.
+ * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
+ * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
+ *
+ * \param [out] thisTupleId the tuple id in \a this corresponding to the returned minimal distance
+ * \param [out] otherTupleId the tuple id in \a other corresponding to the returned minimal distance
+ * \return the minimal distance between the two set of points \a this and \a other.
+ * \sa DataArrayDouble::findClosestTupleId
+ */
+double DataArrayDouble::minimalDistanceTo(const DataArrayDouble *other, int& thisTupleId, int& otherTupleId) const throw(INTERP_KERNEL::Exception)
+{
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> part1=findClosestTupleId(other);
+ int nbOfCompo(getNumberOfComponents());
+ int otherNbTuples(other->getNumberOfTuples());
+ const double *thisPt(begin()),*otherPt(other->begin());
+ const int *part1Pt(part1->begin());
+ double ret=std::numeric_limits<double>::max();
+ for(int i=0;i<otherNbTuples;i++,part1Pt++,otherPt+=nbOfCompo)
+ {
+ double tmp(0.);
+ for(int j=0;j<nbOfCompo;j++)
+ tmp+=(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j])*(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j]);
+ if(tmp<ret)
+ { ret=tmp; thisTupleId=*part1Pt; otherTupleId=i; }
+ }
+ return sqrt(ret);
+}
+
/*!
* This methods returns for each tuple in \a other which tuple in \a this is the closest.
- * So \a this and \a other have to have the same number of components.
+ * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
+ * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
*
* \return a newly allocated (new object to be dealt by the caller) DataArrayInt having \c other->getNumberOfTuples() tuples and one components.
+ * \sa DataArrayDouble::minimalDistanceTo
*/
DataArrayInt *DataArrayDouble::findClosestTupleId(const DataArrayDouble *other) const throw(INTERP_KERNEL::Exception)
{
{
double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2])),zDelta(fabs(bounds[5]-bounds[4]));
double delta=std::max(xDelta,yDelta); delta=std::max(delta,zDelta);
- double characSize=pow((delta*delta*delta)/thisNbOfTuples,1./3.);
+ double characSize=pow((delta*delta*delta)/((double)thisNbOfTuples),1./3.);
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bbox=computeBBoxPerTuple(characSize*1e-12);
- BBTree<3,int> myTree(bbox->getConstPointer(),0,0,getNumberOfTuples(),characSize*1e-12);
- FindClosestTupleIdAlg<3>(myTree,2.*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
+ BBTreePts<3,int> myTree(bbox->getConstPointer(),0,0,getNumberOfTuples(),characSize*1e-12);
+ FindClosestTupleIdAlg<3>(myTree,3.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
break;
}
case 2:
{
double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2]));
double delta=std::max(xDelta,yDelta);
- double characSize=sqrt((delta*delta)/thisNbOfTuples);
+ double characSize=sqrt(delta/(double)thisNbOfTuples);
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bbox=computeBBoxPerTuple(characSize*1e-12);
- BBTree<2,int> myTree(bbox->getConstPointer(),0,0,getNumberOfTuples(),characSize*1e-12);
- FindClosestTupleIdAlg<2>(myTree,2.*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
+ BBTreePts<2,int> myTree(bbox->getConstPointer(),0,0,getNumberOfTuples(),characSize*1e-12);
+ FindClosestTupleIdAlg<2>(myTree,2.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
break;
}
case 1:
{
double characSize=fabs(bounds[1]-bounds[0])/thisNbOfTuples;
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bbox=computeBBoxPerTuple(characSize*1e-12);
- BBTree<1,int> myTree(bbox->getConstPointer(),0,0,getNumberOfTuples(),characSize*1e-12);
- FindClosestTupleIdAlg<1>(myTree,2.*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
+ BBTreePts<1,int> myTree(bbox->getConstPointer(),0,0,getNumberOfTuples(),characSize*1e-12);
+ FindClosestTupleIdAlg<1>(myTree,1.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
break;
}
default:
#include "MEDCouplingTimeLabel.hxx"
#include "MEDCouplingRefCountObject.hxx"
#include "InterpKernelException.hxx"
+#include "BBTreePts.txx"
#include "BBTree.txx"
#include <string>
MEDCOUPLING_EXPORT DataArrayDouble *keepSelectedComponents(const std::vector<int>& compoIds) const throw(INTERP_KERNEL::Exception);
MEDCOUPLING_EXPORT void meldWith(const DataArrayDouble *other) throw(INTERP_KERNEL::Exception);
MEDCOUPLING_EXPORT void findCommonTuples(double prec, int limitTupleId, DataArrayInt *&comm, DataArrayInt *&commIndex) const throw(INTERP_KERNEL::Exception);
+ MEDCOUPLING_EXPORT double minimalDistanceTo(const DataArrayDouble *other, int& thisTupleId, int& otherTupleId) const throw(INTERP_KERNEL::Exception);
MEDCOUPLING_EXPORT DataArrayDouble *duplicateEachTupleNTimes(int nbTimes) const throw(INTERP_KERNEL::Exception);
MEDCOUPLING_EXPORT DataArrayDouble *getDifferentValues(double prec, int limitTupleId=-1) const throw(INTERP_KERNEL::Exception);
MEDCOUPLING_EXPORT DataArrayInt *findClosestTupleId(const DataArrayDouble *other) const throw(INTERP_KERNEL::Exception);
template<int SPACEDIM>
void findCommonTuplesAlg(const double *bbox, int nbNodes, int limitNodeId, double prec, DataArrayInt *c, DataArrayInt *cI) const;
template<int SPACEDIM>
- static void FindClosestTupleIdAlg(const BBTree<SPACEDIM,int>& myTree, double dist, const double *pos, int nbOfTuples, const double *thisPt, int thisNbOfTuples, int *res);
- template<int SPACEDIM>
- static int ComputeBasicClosestTupleIdAlg(const std::vector<int>& elems, const double *thisPt, const double *zePt);
+ static void FindClosestTupleIdAlg(const BBTreePts<SPACEDIM,int>& myTree, double dist, const double *pos, int nbOfTuples, const double *thisPt, int thisNbOfTuples, int *res);
template<int SPACEDIM>
static void FindTupleIdsNearTuplesAlg(const BBTree<SPACEDIM,int>& myTree, const double *pos, int nbOfTuples, double eps,
DataArrayInt *c, DataArrayInt *cI);
}
}
+ PyObject *minimalDistanceTo(const DataArrayDouble *other) const throw(INTERP_KERNEL::Exception)
+ {
+ int thisTupleId,otherTupleId;
+ double r0=self->minimalDistanceTo(other,thisTupleId,otherTupleId);
+ PyObject *ret=PyTuple_New(3);
+ PyTuple_SetItem(ret,0,PyFloat_FromDouble(r0));
+ PyTuple_SetItem(ret,1,PyInt_FromLong(thisTupleId));
+ PyTuple_SetItem(ret,2,PyInt_FromLong(otherTupleId));
+ return ret;
+ }
+
PyObject *getMaxValue() const throw(INTERP_KERNEL::Exception)
{
int tmp;
