From: ageay <ageay>
Date: Tue, 12 Mar 2013 13:01:42 +0000 (+0000)
Subject: DataArrayDouble::findClosestTupleId (for Gauss<->Gauss interpolation)
X-Git-Tag: V6_main_FINAL~304
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=0e1a76b83915a530648b0baf6262160739f8eea6;p=tools%2Fmedcoupling.git

DataArrayDouble::findClosestTupleId (for Gauss<->Gauss interpolation)
---

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 <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
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<SPACEDIM,int>& myTr
 }
 
 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; }
         }
     }
 }
@@ -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<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)
 {
@@ -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<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:
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 <string>
@@ -213,6 +214,7 @@ namespace ParaMEDMEM
     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);
@@ -317,9 +319,7 @@ namespace ParaMEDMEM
     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);
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;