--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java	Fri Jan 22 18:22:11 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java	Fri Jan 22 18:51:47 2010 +0000
@@ -1,17 +1,16 @@
 package de.intevation.gnv.math;
 
-import java.util.List;
-import java.util.Arrays;
-import java.util.Collections;
-
-import java.io.Serializable;
-
-import java.awt.Dimension;
-
 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.Arrays;
+import java.util.List;
 
 import org.apache.log4j.Logger;
 
@@ -26,6 +25,8 @@
     public static final int DEFAULT_WIDTH  = 1024;
     public static final int DEFAULT_HEIGHT =  768;
 
+    public static final double EPS = 1e-6d;
+
     protected int width;
     protected int height;
 
@@ -94,20 +95,18 @@
             return false;
         }
 
-        double [] buffer = Interpolation2D.calculateBuffer(points);
-        double dxMax = buffer[0];
-        double dyMax = buffer[1];
+        List<GridCell> cells = GridCell.pointsToGridCells(points);
 
-        if (debug) {
-            log.debug("buffer size: " + dxMax + " / " + dyMax);
+        if (cells.isEmpty()) {
+            log.warn("no cells found");
+            return false;
         }
 
         // 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(
@@ -147,8 +146,8 @@
         double [] raster = new double[width*height];
         Arrays.fill(raster, Double.NaN);
 
-        Envelope queryBuffer  = new Envelope();
-        XYColumn [] neighbors = new XYColumn[4];
+        Envelope            queryBuffer  = new Envelope();
+        GridCell.CellFinder finder       = new GridCell.CellFinder();
 
         i = 0;
         for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
@@ -159,72 +158,53 @@
             linearToMap.locate(p, center);
 
             queryBuffer.init(
-                center.x - dxMax, center.x + dxMax, 
-                center.y - dyMax, center.y + dyMax);
-
-            List potential = spatialIndex.query(queryBuffer);
-
-            L1Comparator invL1 = new L1Comparator(center);
-            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 j = potential.size()-1; j >= 0; --j) {
-                XYColumn n = (XYColumn)potential.get(j);
-                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) {
+                continue;
             }
 
-            int numNeighbors = Integer.bitCount(mask);
+            XYColumn n0 = (XYColumn)found.p1;
+            XYColumn n1 = (XYColumn)found.p2;
+            XYColumn n2 = (XYColumn)found.p3;
+            XYColumn n3 = (XYColumn)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])
-            &&  neighbors[0].prepare(xyDepth)
-            &&  neighbors[1].prepare(xyDepth)
-            &&  neighbors[2].prepare(xyDepth)
-            &&  neighbors[3].prepare(xyDepth)
+            if (n0.prepare(xyDepth)
+            &&  n1.prepare(xyDepth)
+            &&  n2.prepare(xyDepth)
+            &&  n3.prepare(xyDepth)
             ) {
                 double y1 = Interpolation2D.interpolate(
-                    neighbors[0].x, neighbors[0].y,
-                    neighbors[1].x, neighbors[1].y,
+                    n0.x, n0.y,
+                    n1.x, n1.y,
                     center.x);
                 double y2 = Interpolation2D.interpolate(
-                    neighbors[2].x, neighbors[2].y,
-                    neighbors[3].x, neighbors[3].y,
+                    n2.x, n2.y,
+                    n3.x, n3.y,
                     center.x);
                 int j = i;
                 for (double z = -cellHeight*0.5; 
                     j < raster.length && z >= depth;
                     z -= cellHeight, j += width) {
 
-                    double v0 = neighbors[0].value(z);
-                    double v1 = neighbors[1].value(z);
-                    double v2 = neighbors[2].value(z);
-                    double v3 = neighbors[3].value(z);
+                    double v0 = n0.value(z);
+                    double v1 = n1.value(z);
+                    double v2 = n2.value(z);
+                    double v3 = n3.value(z);
 
                     double z1 = Interpolation2D.interpolate(
-                        neighbors[0].x, v0,
-                        neighbors[1].x, v1,
+                        n0.x, v0,
+                        n1.x, v1,
                         center.x);
                     double z2 = Interpolation2D.interpolate(
-                        neighbors[2].x, v2,
-                        neighbors[3].x, v3,
+                        n2.x, v2,
+                        n3.x, v3,
                         center.x);
                     double value = Interpolation2D.interpolate(
                         y1, z1,