diff --git a/Geo/closestPoint.cpp b/Geo/closestPoint.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..67d43a41020abc5f179de93cec51073d55de3588
--- /dev/null
+++ b/Geo/closestPoint.cpp
@@ -0,0 +1,69 @@
+#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]);
+}
+
diff --git a/Geo/closestPoint.h b/Geo/closestPoint.h
new file mode 100644
index 0000000000000000000000000000000000000000..24ecc73fa4805abf2a2a437048883f6f86f76ee8
--- /dev/null
+++ b/Geo/closestPoint.h
@@ -0,0 +1,26 @@
+#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
diff --git a/Numeric/hausdorffDistance.cpp b/Numeric/hausdorffDistance.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0c7ba1af3072bf40590bd42fe8d44095edd4c1aa
--- /dev/null
+++ b/Numeric/hausdorffDistance.cpp
@@ -0,0 +1,115 @@
+/*
+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));
+}
diff --git a/Numeric/hausdorffDistance.h b/Numeric/hausdorffDistance.h
new file mode 100644
index 0000000000000000000000000000000000000000..09f4ef39a213d979bf0dd86aa51a5acbc0bd0b88
--- /dev/null
+++ b/Numeric/hausdorffDistance.h
@@ -0,0 +1,7 @@
+#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