--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Wed Jan 20 08:50:10 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Wed Jan 20 14:47:30 2010 +0000
@@ -1,19 +1,24 @@
 package de.intevation.gnv.math;
 
-import java.util.ArrayList;
-import java.util.List;
-import java.util.HashMap;
-import java.util.Collections;
-
 import com.vividsolutions.jts.geom.Coordinate;
 import com.vividsolutions.jts.geom.Envelope;
 
 import com.vividsolutions.jts.index.quadtree.Quadtree;
 
+import gnu.trove.TDoubleArrayList;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.List;
+
+import org.apache.commons.math.stat.descriptive.moment.Mean;
+import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;
+
 import org.apache.log4j.Logger;
 
 /**
- *  @author Sascha L. Teichmann
+ *  @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
  */
 public final class Interpolation2D
 {
@@ -26,6 +31,108 @@
     private Interpolation2D() {
     }
 
+    private static final double dist(double a, double b) {
+        a -= b;
+        return Math.sqrt(a*a);
+    }
+
+    private static final double OUTLIER_EPS = 0.4d;
+
+    public static final double [] calculateBufferRemoveOutlier(
+        List <? extends Point2d> points
+    ) {
+        HashMap<Integer, ArrayList<Point2d>> iMap = 
+            new HashMap<Integer, ArrayList<Point2d>>();
+
+        HashMap<Integer, ArrayList<Point2d>> jMap = 
+            new HashMap<Integer, ArrayList<Point2d>>();
+
+        for (int k = points.size()-1; k >= 0; --k) {
+            Point2d p = points.get(k);
+
+            ArrayList<Point2d> jList = jMap.get(p.j);
+            ArrayList<Point2d> iList = iMap.get(p.i);
+
+            if (jList == null) {
+                jMap.put(p.j, jList = new ArrayList<Point2d>());
+            }
+            jList.add(p);
+
+            if (iList == null) {
+                iMap.put(p.i, iList = new ArrayList<Point2d>());
+            }
+            iList.add(p);
+        }
+
+        TDoubleArrayList deltas = new TDoubleArrayList();
+
+        Mean mean            = new Mean();
+        StandardDeviation sd = new StandardDeviation();
+
+        for (ArrayList<Point2d> v: jMap.values()) {
+            Collections.sort(v, Point2d.Y_COMPARATOR);
+            for (int i = 1, L = v.size(); i < L; ++i) {
+                double dy = Math.abs(v.get(i).y - v.get(i-1).y);
+                deltas.add(dy);
+                mean.increment(dy);
+                sd  .increment(dy);
+            }
+        }
+
+        deltas.sort();
+
+        double m   = mean.getResult();
+        double s   = sd  .getResult();
+        double eps = OUTLIER_EPS*s;
+
+        int yi = deltas.size() - 1;
+        while (yi > 0 && dist(deltas.get(yi), m) > eps) {
+            --yi;
+        }
+
+        double dyMax = deltas.get(yi) + 1e-5d;
+
+        if (log.isDebugEnabled()) {
+            log.debug("mean  y: " + m);
+            log.debug("sigma y: " + s);
+        }
+
+        deltas.reset();
+        mean.clear();
+        sd  .clear();
+
+        for (ArrayList<Point2d> v: iMap.values()) {
+            Collections.sort(v, Point2d.X_COMPARATOR);
+            for (int i = 1, L = v.size(); i < L; ++i) {
+                double dx = Math.abs(v.get(i).x - v.get(i-1).x);
+                deltas.add(dx);
+                mean.increment(dx);
+                sd  .increment(dx);
+            }
+        }
+
+        deltas.sort();
+
+        m   = mean.getResult();
+        s   = sd  .getResult();
+        eps = OUTLIER_EPS*s;
+
+        int xi = deltas.size() - 1;
+
+        while (xi > 0 && dist(deltas.get(xi), m) > eps) {
+            --xi;
+        }
+
+        double dxMax = deltas.get(xi) + 1e-5d;
+
+        if (log.isDebugEnabled()) {
+            log.debug("mean  y: " + m);
+            log.debug("sigma y: " + s);
+        }
+
+        return new double [] { dxMax, dyMax };
+    }
+
     public static final double [] calculateBuffer(
         List <? extends Point2d> points
     ) {