--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java	Fri Jan 22 15:45:59 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java	Fri Jan 22 18:22:11 2010 +0000
@@ -3,15 +3,13 @@
 import com.vividsolutions.jts.geom.Coordinate;
 import com.vividsolutions.jts.geom.Envelope;
 
-import com.vividsolutions.jts.index.quadtree.Quadtree;
+import com.vividsolutions.jts.index.strtree.STRtree;
 
 import java.awt.Dimension;
 
 import java.io.Serializable;
 
-import java.util.ArrayList;
 import java.util.Arrays;
-import java.util.Collections;
 import java.util.List;
 
 import org.apache.log4j.Logger;
@@ -44,6 +42,9 @@
         return raster;
     }
 
+    public static final double EPS = 1e-6d;
+
+
     public boolean interpolate(
         List<? extends Point2d> points,
         Envelope                boundingBox,
@@ -53,9 +54,14 @@
         boolean debug = log.isDebugEnabled();
 
         if (points == null || points.isEmpty()) {
-            if (debug) {
-                log.debug("no points to interpolate");
-            }
+            log.warn("no points to interpolate");
+            return false;
+        }
+
+        List<GridCell> cells = GridCell.pointsToGridCells(points);
+
+        if (cells.isEmpty()) {
+            log.warn("no cells to interpolate");
             return false;
         }
 
@@ -64,47 +70,25 @@
 
         double cellWidth  = boundingBox.getWidth()  / W;
         double cellHeight = boundingBox.getHeight() / H;
-        Envelope gridEnvelope = new Envelope(boundingBox);
-
-        gridEnvelope.expandBy(2d*cellWidth, 2d*cellWidth);
-
-        ArrayList<Point2d>relevantPoints =
-            new ArrayList<Point2d>();
 
-        for (int i = points.size()-1; i >= 0; --i) {
-            if (gridEnvelope.contains(points.get(i))) {
-                relevantPoints.add(points.get(i));
-            }
-        }
+        STRtree spatialIndex = new STRtree();
 
-        double [] buffer = Interpolation2D
-            .calculateBufferRemoveOutlier(relevantPoints);
-        double dxMax = buffer[0];
-        double dyMax = buffer[1];
-
-        Quadtree spatialIndex = new Quadtree();
-
-        for (int i = relevantPoints.size()-1; i >= 0; --i) {
-            Point2d p = relevantPoints.get(i);
-            spatialIndex.insert(p.envelope(), p);
+        for (GridCell cell: cells) {
+            spatialIndex.insert(cell.getEnvelope(), cell);
         }
 
         if (debug) {
-            log.debug("points: "        + relevantPoints.size());
             log.debug("width:  "        + boundingBox.getWidth());
             log.debug("height:  "       + boundingBox.getHeight());
             log.debug("sample width:  " + W);
             log.debug("sample height: " + H);
             log.debug("cell width:  "   + cellWidth);
             log.debug("cell height: "   + cellHeight);
-            log.debug("buffer x: "      + dxMax);
-            log.debug("buffer y: "      + dyMax);
         }
 
-        Envelope     queryBuffer = new Envelope();
-        Point2d []   neighbors   = new Point2d[4];
-        Coordinate   center      = new Coordinate();
-        L1Comparator invL1       = new L1Comparator(center);
+        Envelope   queryBuffer     = new Envelope();
+        Coordinate center          = new Coordinate();
+        GridCell.CellFinder finder = new GridCell.CellFinder();
 
         double [] raster = new double[W*H];
         Arrays.fill(raster, Double.NaN);
@@ -119,64 +103,38 @@
             double y = j*cellHeight + 0.5d*cellHeight + minY;
             double x = 0.5d*cellWidth + minX;
             for (int i = 0; i < W; ++i, x += cellWidth) {
-                queryBuffer.init(
-                    x - dxMax, x + dxMax, 
-                    y - dyMax, y + dyMax);
 
-                List potential = spatialIndex.query(queryBuffer);
+                queryBuffer.init(x - EPS, x + EPS, y - EPS, y + EPS);
+                center.x = x; center.y = y;
+                finder.prepare(center);
+                spatialIndex.query(queryBuffer, finder);
 
-                if (potential.size() < 4) {
+                GridCell found = finder.found;
+
+                if (found == null) {
                     continue;
                 }
 
-                center.x = x; center.y = y;
-                Collections.sort(potential, invL1);
-
-                neighbors[0] = neighbors[1] = 
-                neighbors[2] = neighbors[3] = null;
-
-                int mask = 0;
-                // reversed order is essential here!
-                for (int k = potential.size()-1; k >= 0; --k) {
-                    Point2d n = (Point2d)potential.get(k);
-                    int code = n.x > center.x ? 1 : 0;
-                    if (n.y > center.y) code |= 2;
-                    neighbors[code] = n;
-                    mask |= 1 << code;
-                }
-
-                int numNeighbors = Integer.bitCount(mask);
-
-                // only interpolate if we have all four neighbors
-                // we do not have any gaps and we are below surface.
-                if (numNeighbors == 4
-                && !neighbors[0].hasIGap(neighbors[1])
-                && !neighbors[1].hasJGap(neighbors[3])
-                && !neighbors[3].hasIGap(neighbors[2])
-                && !neighbors[2].hasJGap(neighbors[0])
-                && depth.depth(center) <= 0d
-                ) {
-                    double z1 = Interpolation2D.interpolate(
-                        neighbors[0].x, neighbors[0].z,
-                        neighbors[1].x, neighbors[1].z,
-                        center.x);
-                    double z2 = Interpolation2D.interpolate(
-                        neighbors[2].x, neighbors[2].z,
-                        neighbors[3].x, neighbors[3].z,
-                        center.x);
-                    double y1 = Interpolation2D.interpolate(
-                        neighbors[0].x, neighbors[0].y,
-                        neighbors[1].x, neighbors[1].y,
-                        center.x);
-                    double y2 = Interpolation2D.interpolate(
-                        neighbors[2].x, neighbors[2].y,
-                        neighbors[3].x, neighbors[3].y,
-                        center.x);
-                    raster[row+i] = Interpolation2D.interpolate(
-                        y1, z1,
-                        y2, z2,
-                        center.y);
-                }
+                double z1 = Interpolation2D.interpolate(
+                    found.p1.x, found.p1.z,
+                    found.p2.x, found.p2.z,
+                    center.x);
+                double z2 = Interpolation2D.interpolate(
+                    found.p3.x, found.p3.z,
+                    found.p4.x, found.p4.z,
+                    center.x);
+                double y1 = Interpolation2D.interpolate(
+                    found.p1.x, found.p1.y,
+                    found.p2.x, found.p2.y,
+                    center.x);
+                double y2 = Interpolation2D.interpolate(
+                    found.p3.x, found.p3.y,
+                    found.p4.x, found.p4.y,
+                    center.x);
+                raster[row+i] = Interpolation2D.interpolate(
+                    y1, z1,
+                    y2, z2,
+                    center.y);
             }
         }