From 0e1a76b83915a530648b0baf6262160739f8eea6 Mon Sep 17 00:00:00 2001 From: ageay Date: Tue, 12 Mar 2013 13:01:42 +0000 Subject: [PATCH] DataArrayDouble::findClosestTupleId (for Gauss<->Gauss interpolation) --- src/INTERP_KERNEL/BBTree.txx | 6 +- src/INTERP_KERNEL/BBTreePts.txx | 221 +++++++++++++++++++++++ src/INTERP_KERNEL/Makefile.am | 1 + src/MEDCoupling/MEDCouplingMemArray.cxx | 89 +++++---- src/MEDCoupling/MEDCouplingMemArray.hxx | 6 +- src/MEDCoupling_Swig/MEDCouplingCommon.i | 11 ++ 6 files changed, 292 insertions(+), 42 deletions(-) create mode 100644 src/INTERP_KERNEL/BBTreePts.txx diff --git a/src/INTERP_KERNEL/BBTree.txx b/src/INTERP_KERNEL/BBTree.txx index 7b3d03782..dfb54606d 100644 --- a/src/INTERP_KERNEL/BBTree.txx +++ b/src/INTERP_KERNEL/BBTree.txx @@ -52,7 +52,9 @@ public: \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 @@ -64,7 +66,7 @@ public: \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) diff --git a/src/INTERP_KERNEL/BBTreePts.txx b/src/INTERP_KERNEL/BBTreePts.txx new file mode 100644 index 000000000..934666e2f --- /dev/null +++ b/src/INTERP_KERNEL/BBTreePts.txx @@ -0,0 +1,221 @@ +// 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 +#include + +#include +#include +#include + +template +class BBTreePts +{ + +private: + BBTreePts* _left; + BBTreePts* _right; + int _level; + double _max_left; + double _min_right; + const double *_pts; + typename std::vector _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(nodes, nodes+nbelems/2, nodes+nbelems); + double median=*(nodes+nbelems/2); + delete [] nodes; + std::vector 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::max(); + double min_right= std::numeric_limits::max(); + for(int i=0;imedian) + { + new_elems_right.push_back(elem); + if(mxmax_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::max(); + for(ConnType i=0;i<_nbelems;i++) + { + const double* const bb_ptr=_pts+_elems[i]*dim; + double tmp=0.; + for(int idim=0;idimgetElementsAroundPoint2(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& 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_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 diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index a21327824..bb81f0a65 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -28,6 +28,7 @@ lib_LTLIBRARIES = libinterpkernel.la salomeinclude_HEADERS = \ BBTree.txx \ +BBTreePts.txx \ BoundingBox.hxx \ CellModel.hxx \ ConvexIntersector.hxx \ diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index 04e3cff2a..69ca3dd43 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -85,42 +85,26 @@ void DataArrayDouble::FindTupleIdsNearTuplesAlg(const BBTree& myTr } template -int DataArrayDouble::ComputeBasicClosestTupleIdAlg(const std::vector& elems, const double *thisPt, const double *zePt) -{ - double val=std::numeric_limits::max(); - int ret=-1; - for(std::vector::const_iterator it=elems.begin();it!=elems.end();it++) - { - double tmp=0.; - for(int j=0;j -void DataArrayDouble::FindClosestTupleIdAlg(const BBTree& myTree, double dist, const double *pos, int nbOfTuples, const double *thisPt, int thisNbOfTuples, int *res) +void DataArrayDouble::FindClosestTupleIdAlg(const BBTreePts& 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::max()); int nbOfTurn=15; - while(nbOfTurn>0) - { - for(int j=0;j elems; - myTree.getIntersectingElems(bbox,elems); nbOfTurn--; - if(elems.empty()) - { failVal=distOpt; distOpt=okVal==std::numeric_limits::max()?2*distOpt:(okVal+failVal)/2.; continue; } - if( elems.size()<=15 || nbOfTurn>0) - { *r=ComputeBasicClosestTupleIdAlg(elems,thisPt,p); break; } + while(true) + { + int elem=-1; + double ret=myTree.getElementsAroundPoint2(p,distOpt,elem); nbOfTurn--; + if(ret!=std::numeric_limits::max()) + { + distOpt=ret>1e-4?ret:dist; + *r=elem; + break; + } else - { okVal=distOpt; distOpt=(distOpt+failVal)/2.; continue; } + { distOpt=2*distOpt; continue; } } } } @@ -1352,11 +1336,42 @@ DataArrayDouble *DataArrayDouble::duplicateEachTupleNTimes(int nbTimes) const th 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 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::max(); + for(int i=0;igetNumberOfTuples() tuples and one components. + * \sa DataArrayDouble::minimalDistanceTo */ DataArrayInt *DataArrayDouble::findClosestTupleId(const DataArrayDouble *other) const throw(INTERP_KERNEL::Exception) { @@ -1381,28 +1396,28 @@ DataArrayInt *DataArrayDouble::findClosestTupleId(const DataArrayDouble *other) { 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 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 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 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: diff --git a/src/MEDCoupling/MEDCouplingMemArray.hxx b/src/MEDCoupling/MEDCouplingMemArray.hxx index b7205f15f..08074596b 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.hxx +++ b/src/MEDCoupling/MEDCouplingMemArray.hxx @@ -25,6 +25,7 @@ #include "MEDCouplingTimeLabel.hxx" #include "MEDCouplingRefCountObject.hxx" #include "InterpKernelException.hxx" +#include "BBTreePts.txx" #include "BBTree.txx" #include @@ -213,6 +214,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT DataArrayDouble *keepSelectedComponents(const std::vector& 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); @@ -317,9 +319,7 @@ namespace ParaMEDMEM template void findCommonTuplesAlg(const double *bbox, int nbNodes, int limitNodeId, double prec, DataArrayInt *c, DataArrayInt *cI) const; template - static void FindClosestTupleIdAlg(const BBTree& myTree, double dist, const double *pos, int nbOfTuples, const double *thisPt, int thisNbOfTuples, int *res); - template - static int ComputeBasicClosestTupleIdAlg(const std::vector& elems, const double *thisPt, const double *zePt); + static void FindClosestTupleIdAlg(const BBTreePts& myTree, double dist, const double *pos, int nbOfTuples, const double *thisPt, int thisNbOfTuples, int *res); template static void FindTupleIdsNearTuplesAlg(const BBTree& myTree, const double *pos, int nbOfTuples, double eps, DataArrayInt *c, DataArrayInt *cI); diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 72ddf3688..144b50669 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -3257,6 +3257,17 @@ namespace ParaMEDMEM } } + 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; -- 2.39.2