diff --git a/Numeric/discreteFrechetDistance.cpp b/Numeric/discreteFrechetDistance.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..10281d41bfa691c9fb3e49eeb9b00e0d953fd1da
--- /dev/null
+++ b/Numeric/discreteFrechetDistance.cpp
@@ -0,0 +1,48 @@
+#include "discreteFrechetDistance.h"
+#include "fullMatrix.h"
+
+static double distance (const SPoint3 &p1, const SPoint3 &p2){
+  return p1.distance(p2);
+} 
+
+static double c(int i, int j,
+		fullMatrix<double> &CA,
+		const std::vector<SPoint3> &P, 
+		const std::vector<SPoint3> &Q) 
+{
+  double CAij;
+  if (CA(i,j)>-1){
+    CAij =  CA(i,j);
+  }
+  else if (i==0 && j==0){
+    CA(i,j) = distance(P[0],Q[0]);     // update the CA permanent
+    CAij = CA(i,j);                    // set the current relevant value
+  }
+  else if (i>0 && j==0){
+    CA(i,j) = std::max(c(i-1,0,CA,P,Q), distance(P[i],Q[1]) );
+    CAij = CA(i,j);
+  }
+  else if (i==0 && j>0){
+    CA(i,j) = std::max( c(0,j-1,CA,P,Q), distance(P[0],Q[j]) );
+    CAij = CA(i,j);
+  }
+  else if (i>0 && j>0){
+    double temp = std::min(std::min(c(i-1,j,CA,P,Q),c(i-1,j-1,CA,P,Q)),c(i,j-1,CA,P,Q));
+    CA(i,j) = std::max(temp,distance(P[i],Q[j]));
+    CAij = CA(i,j);
+  }
+  else{
+    CA(i,j) = 1.e22;
+  }
+  return CAij;
+}
+
+double discreteFrechetDistance (const std::vector<SPoint3> &P, 
+				const std::vector<SPoint3> &Q){
+  const int sP = P.size();
+  const int sQ = Q.size();
+  fullMatrix<double> CA (sP,sQ);
+  CA.setAll(-1.0);
+  double cm = c(sP-1,sQ-1,CA,P,Q);
+  return cm;
+}
diff --git a/Numeric/discreteFrechetDistance.h b/Numeric/discreteFrechetDistance.h
new file mode 100644
index 0000000000000000000000000000000000000000..439b0c9b0ed34daa2377161c4574b46a49689ab6
--- /dev/null
+++ b/Numeric/discreteFrechetDistance.h
@@ -0,0 +1,7 @@
+#ifndef _DISCRETE_FRECHET_DISTANCE_
+#define _DISCRETE_FRECHET_DISTANCE_
+#include <vector>
+#include "SPoint3.h"
+double discreteFrechetDistance (const std::vector<SPoint3> &P, 
+				const std::vector<SPoint3> &Q);
+#endif