oops

Geo/closestPoint.cpp

+#include "closestPoint.h"
+#include "GEntity.h"
+#include "GEdge.h"
+#include "GFace.h"
+#include <vector>
+
+static void oversample (std::vector<SPoint3> &s, double tol){
+  std::vector<SPoint3> t;
+  for (unsigned int i=1;i<s.size();i++){
+    SPoint3 p0 = s[i-1];
+    SPoint3 p1 = s[i];
+    double d = p0.distance(p1);
+    int N = (int) (d / tol);
+    t.push_back(p0);
+    for (int j=1;j<N;j++){
+      const double xi = (double) j/ N;
+      t.push_back(p0 + (p1-p0)*xi);
+    }
+    t.push_back(p1);
+  }
+  s = t;
+}
+
+closestPointFinder :: closestPointFinder (GEntity *ge, double e) : _ge (ge) , _tolerance (e) 
+{
+#if defined(HAVE_ANN)
+  std::vector<SPoint3> pts;
+  if (ge->dim() == 1){
+    GEdge *edge = ge->cast2Edge();
+    if (!edge)Msg::Fatal("in the constructor of closestPoint");
+    std::vector<double> ts;
+    edge->discretize(_tolerance, pts,ts);
+    //    printf("%d points\n",pts.size());
+    oversample (pts, _tolerance);
+    //    printf("%d points after oversampling\n",pts.size());
+  }
+  index = new ANNidx[1];
+  dist = new ANNdist[1];
+  zeronodes = annAllocPts(pts.size(), 3);
+  for (unsigned int k=0;k<pts.size();k++){
+    zeronodes[k][0] = pts[k].x();
+    zeronodes[k][1] = pts[k].y();
+    zeronodes[k][2] = pts[k].z();
+  }
+  kdtree = new ANNkd_tree(zeronodes, pts.size(), 3);
+#else
+  Msg::Fatal("Gmsh should be compiled using ANN");
+#endif
+}
+
+closestPointFinder :: ~closestPointFinder ()
+{
+#if defined(HAVE_ANN)
+  if(kdtree) delete kdtree;
+  if(zeronodes) annDeallocPts(zeronodes);
+  delete[]index;
+  delete[]dist;  
+#endif
+}
+
+SPoint3 closestPointFinder ::operator() (const SPoint3 &p)
+{
+  double xyz[3] = {p.x(),p.y(),p.z()};
+  kdtree->annkSearch(xyz, 1, index, dist);
+  return SPoint3(zeronodes[index[0]][0],
+		 zeronodes[index[0]][1],
+		 zeronodes[index[0]][2]);  
+}
+

Geo/closestPoint.h

+#ifndef _CLOSEST_POINT_H_
+#define _CLOSEST_POINT_H_
+#include "GmshConfig.h"
+#if defined(HAVE_ANN)
+#include "ANN/ANN.h"
+#endif
+#include "SPoint3.h"
+class GEntity;
+class closestPointFinder
+{
+#if defined(HAVE_ANN)
+  ANNkd_tree *kdtree;
+  ANNpointArray zeronodes;
+  ANNidxArray index;
+  ANNdistArray dist;
+#endif
+  GEntity* _ge;
+  double _tolerance;
+ public :
+  closestPointFinder (GEntity*, double);
+  ~closestPointFinder ();
+  SPoint3 operator () (const SPoint3 & p);
+  inline double tol() const {return _tolerance;}
+};
+
+#endif

Numeric/hausdorffDistance.cpp

+/*
+compute the hausdorff distance between two polygonal curves 
+in n*m time where n and m are the nb of points of the 
+polygonal curves
+*/
+#include "SVector3.h"
+#include "hausdorffDistance.h"
+
+static double intersect (SPoint3 &q, SVector3 &n, // a plane
+			 SPoint3 &p1, SPoint3 &p2, // a segment
+			 SPoint3 &result)
+{
+  // n * (x-q) = 0
+  // x = p1 + t (p2-p1)
+  // n *(p1 + t (p2-p1) - q) = 0
+  // t = n*(q-p1) / n*(p2-p1)
+
+  const double t = dot(n,q-p1)/dot(n,p2-p1);
+  result = p1 *(1.-t)+ p2*t;
+  return t;
+}
+
+static double projection (SPoint3 &p1, SPoint3 &p2, SPoint3 &q, SPoint3 &result){
+  // x = p1 + t (p2 - p1)
+  // (x - q) * (p2 - p1) = 0
+  // (p1 + t (p2 - p1) - q) (p2 - p1) = 0
+  // t = (q-p1) * (p2-p1) / (p2-p1)^2
+  SVector3 p21 = p2-p1;
+  const double t = dot(q-p1,p21)/dot(p21,p21);
+  result = p1 *(1.-t)+ p2*t;
+  return t;
+}  
+
+static SPoint3 closestPoint (SPoint3 &p1, SPoint3 &p2, SPoint3 &q){
+  double result;
+  const double t = projection (p1,p2,q,result);
+  if (t >= 0.0 && t <= 1.0)return result;
+  if (t < 0)return p1;
+  return p2;
+}  
+
+
+double closestPoint (const std::vector<SPoint3> &P, const SPoint3 &p, SPoint3 & result){
+  double closestDistance = 1.e22;
+  for (unsigned int i=1;i<P.size();i++){
+    SPoint3 q = closestPoint (P[i-1],P[i],p);
+    const double pq = p.distance(q);
+    if (pq < closestDistance){
+      closestDistance = pq;
+      result = q;
+    }   
+  }
+  return closestDistance;
+} 
+
+// we test all points of P plus all points that are the intersections
+// of angle bissectors of Q with P
+double oneSidedHausdorffDistance (const std::vector<SPoint3> &P, 
+				  const std::vector<SPoint3> &Q, 
+				  SPoint3 &p1, SPoint3 &p2){
+  const double hausdorffDistance = 0.0;
+  
+  // first test the points
+  for (unsigned int i=0;i<P.size();i++){
+    SPoint3 result;
+    double d = closestPoint (Q, P[i], result);
+    if (d > hausdorffDistance){
+      hausdorffDistance = d;
+      p1 = P[i];
+      p2 = result;
+    }
+  }
+  // compute angle bissectors intersections
+  std::vector<SPoint3> intersections;
+  for (unsigned int i=1;i<Q.size()-1;i++){
+    SPoint3 a = Q[i-1];
+    SPoint3 b = Q[i];
+    SPoint3 c = Q[i+1];
+    SVector3 ba = b-a;
+    SVector3 ca = c-a;
+    SVector3 bissector = (ba+ca);
+    SVector3 n; // normal to the bissector plane
+    if (bissector.norm == 0){
+      ba.normalize();
+      n = ba;
+    }
+    else{
+      SVector3 b = crossprod (bissector,ba);
+      n = crossprod (b,bissector);
+      n.normalize();
+    }
+    for (unsigned int i=1;i<P.size();i++){
+      SPoint3 result;
+      const double t = intersect (b, n, P[i-1],P[i],result);
+      if (t >=0 && t <=1)intersections.push_back(result);	
+    }
+  }
+
+  for (unsigned int i=0;i<intersections.size();i++){
+    SPoint3 result;
+    double d = closestPoint (Q, intersections[i], result);
+    if (d > hausdorffDistance){
+      hausdorffDistance = d;
+      p1 = P[i];
+      p2 = result;
+    }
+  }
+  return hausdorffDistance;
+}
+
+double hausdorffDistance (const std::vector<SPoint3> &P, 
+			  const std::vector<SPoint3> &Q){
+  return std::max(oneSidedHausdorffDistance (P,Q),
+		  oneSidedHausdorffDistance (Q,P));  
+}

Numeric/hausdorffDistance.h

+#ifndef _HAUSDORFF_DISTANCE_H_
+#define _HAUSDORFF_DISTANCE_H_
+#include <vector>
+#include "SPoint3.h"
+double hausdorffDistance (const std::vector<SPoint3> &P, 
+			  const std::vector<SPoint3> &Q);
+#endif