--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Sat Jan 23 18:10:34 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Sat Jan 23 21:16:45 2010 +0000
@@ -3,18 +3,10 @@
 import com.vividsolutions.jts.geom.Coordinate;
 import com.vividsolutions.jts.geom.Envelope;
 
-import com.vividsolutions.jts.index.quadtree.Quadtree;
-
-import gnu.trove.TDoubleArrayList;
+import com.vividsolutions.jts.index.strtree.STRtree;
 
-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;
 
 /**
@@ -24,6 +16,11 @@
 {
     private static Logger log = Logger.getLogger(Interpolation2D.class);
 
+    public static final int CULL_POINT_THRESHOLD = Integer.getInteger(
+        "gnv.interpolation2d.cull.point.threshold", 1000);
+
+    public static final double EPS = 1e-6d;
+
     public interface Consumer {
         void interpolated(Coordinate point, boolean success);
     } // interface Consumer
@@ -31,164 +28,22 @@
     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
+    public static Envelope pathBoundingBox(
+        List<? extends Coordinate> path,
+        double extra
     ) {
-        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);
+        int N = path.size();
+        Envelope area = new Envelope(path.get(N-1));
 
-            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);
-                double dy = v.get(i).distance(v.get(i-1));
-                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);
+        for (int i = N-2; i >= 0; --i) {
+            area.expandToInclude(path.get(i));
         }
 
-        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);
-                double dx = v.get(i).distance(v.get(i-1));
-                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
-    ) {
-        HashMap<Integer, ArrayList<Point2d>> iMap = 
-            new HashMap<Integer, ArrayList<Point2d>>();
+        area.expandBy(
+            extra*area.getWidth(),
+            extra*area.getHeight());
 
-        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);
-        }
-
-        double dxMax = -Double.MAX_VALUE;
-        double dyMax = -Double.MAX_VALUE;
-
-        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);
-                if (dy > dyMax) {
-                    dyMax = dy;
-                }
-            }
-        }
-
-        dyMax += 1e-5d;
-
-        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);
-                if (dx > dxMax) {
-                    dxMax = dx;
-                }
-            }
-        }
-
-        dxMax += 1e-5d;
-
-        return new double [] { dxMax, dyMax };
+        return area;
     }
 
     public static void interpolate(
@@ -214,20 +69,23 @@
             return;
         }
 
-        double [] buffer = calculateBuffer(points);
-        double dxMax = buffer[0];
-        double dyMax = buffer[1];
+        Envelope relevantArea = M > CULL_POINT_THRESHOLD
+			? pathBoundingBox(path, 0.05d)
+			: null;
+        
+        List<GridCell> cells = GridCell.pointsToGridCells(
+            points, relevantArea);
 
-        if (debug) {
-            log.debug("buffer size: " + dxMax + " / " + dyMax);
+        if (cells.isEmpty()) {
+            log.warn("no cells found");
+            return;
         }
 
         // put into spatial index to speed up finding neighbors.
-        Quadtree spatialIndex = new Quadtree();
+        STRtree spatialIndex = new STRtree();
 
-        for (int i = 0; i < M; ++i) {
-            Point2d p = points.get(i);
-            spatialIndex.insert(p.envelope(), p);
+        for (GridCell cell: cells) {
+            spatialIndex.insert(cell.getEnvelope(), cell);
         }
 
         LinearToMap linearToMap = new LinearToMap(
@@ -235,85 +93,59 @@
 
         double dP = (to - from)/steps;
 
-        Coordinate center = new Coordinate();
-
-        Envelope queryBuffer = new Envelope();
-
-        Point2d [] neighbors = new Point2d[4];
+        Coordinate          center      = new Coordinate();
+        Envelope            queryBuffer = new Envelope();
+        GridCell.CellFinder finder      = new GridCell.CellFinder();
 
         int missedInterpolations = 0;
-        int interpolations = 0;
-
-        L1Comparator invL1 = new L1Comparator(center);
+        int interpolations       = 0;
 
         for (double p = from; p <= to; p += dP) {
             if (!linearToMap.locate(p, center)) {
                 continue;
             }
             queryBuffer.init(
-                center.x - dxMax, center.x + dxMax, 
-                center.y - dyMax, center.y + dyMax);
-
-            List potential = spatialIndex.query(queryBuffer);
-
-            Collections.sort(potential, invL1);
-
-            neighbors[0] = neighbors[1] = 
-            neighbors[2] = neighbors[3] = null;
+                center.x - EPS, center.x + EPS, 
+                center.y - EPS, center.y + EPS);
 
-            /* bit code of neighbors
-               0---1
-               | x |
-               2---3
-            */
+            finder.prepare(center);
+            spatialIndex.query(queryBuffer, finder);
 
-            int mask = 0;
-            // reversed order is essential here!
-            for (int i = potential.size()-1; i >= 0; --i) {
-                Point2d n = (Point2d)potential.get(i);
-                int code = n.x > center.x ? 1 : 0;
-                if (n.y > center.y) code |= 2;
-                neighbors[code] = n;
-                mask |= 1 << code;
+            GridCell found = finder.found;
+
+            if (found == null) {
+                consumer.interpolated(center, false);
+                ++missedInterpolations;
+                continue;
             }
 
-            int numNeighbors = Integer.bitCount(mask);
+            Point2d n0 = found.p1;
+            Point2d n1 = found.p2;
+            Point2d n2 = found.p3;
+            Point2d n3 = found.p4;
 
-            // only interpolate if we have all four neighbors
-            // and we do not have any gaps.
-            if (numNeighbors == 4
-            && !neighbors[0].hasIGap(neighbors[1])
-            && !neighbors[1].hasJGap(neighbors[3])
-            && !neighbors[3].hasIGap(neighbors[2])
-            && !neighbors[2].hasJGap(neighbors[0])
-            ) {
-                double z1 = interpolate(
-                    neighbors[0].x, neighbors[0].z,
-                    neighbors[1].x, neighbors[1].z,
-                    center.x);
-                double z2 = interpolate(
-                    neighbors[2].x, neighbors[2].z,
-                    neighbors[3].x, neighbors[3].z,
-                    center.x);
-                double y1 = interpolate(
-                    neighbors[0].x, neighbors[0].y,
-                    neighbors[1].x, neighbors[1].y,
-                    center.x);
-                double y2 = interpolate(
-                    neighbors[2].x, neighbors[2].y,
-                    neighbors[3].x, neighbors[3].y,
-                    center.x);
-                center.z = interpolate(
-                    y1, z1,
-                    y2, z2,
-                    center.y);
-                consumer.interpolated(center, true);
-                ++interpolations;
-            }
-            else {
-                consumer.interpolated(center, false);
-                ++missedInterpolations;
-            }
+            double z1 = interpolate(
+                n0.x, n0.z,
+                n1.x, n1.z,
+                center.x);
+            double z2 = interpolate(
+                n2.x, n2.z,
+                n3.x, n3.z,
+                center.x);
+            double y1 = interpolate(
+                n0.x, n0.y,
+                n1.x, n1.y,
+                center.x);
+            double y2 = interpolate(
+                n2.x, n2.y,
+                n3.x, n3.y,
+                center.x);
+            center.z = interpolate(
+                y1, z1,
+                y2, z2,
+                center.y);
+            consumer.interpolated(center, true);
+            ++interpolations;
         }
 
         if (debug) {