--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java	Thu Jan 14 10:34:05 2010 +0000
@@ -0,0 +1,156 @@
+package de.intevation.gnv.math;
+
+import java.io.Serializable;
+
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.List;
+import java.util.HashMap;
+import java.util.Collections;
+
+import java.awt.Dimension;
+
+import com.vividsolutions.jts.geom.Coordinate;
+import com.vividsolutions.jts.geom.Envelope;
+
+import com.vividsolutions.jts.index.quadtree.Quadtree;
+
+import org.apache.log4j.Logger;
+
+/**
+ *  @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
+ */
+public class AreaInterpolation
+implements   Serializable
+{
+    private static Logger log = Logger.getLogger(AreaInterpolation.class);
+
+    protected double [] raster;
+
+    protected int width;
+    protected int height;
+
+    public AreaInterpolation() {
+    }
+
+    public int getWidth() {
+        return width;
+    }
+
+    public int getHeight() {
+        return height;
+    }
+
+    public double [] getRaster() {
+        return raster;
+    }
+
+    public boolean interpolate(
+        List<? extends Point2d> points,
+        Envelope                boundingBox,
+        Dimension               samples,
+        XYDepth                 depth
+    ) {
+        if (points == null || points.isEmpty()) {
+            return false;
+        }
+
+        double [] buffer = Interpolation2D.calculateBuffer(points);
+        double dxMax = buffer[0];
+        double dyMax = buffer[1];
+
+        Quadtree spatialIndex = new Quadtree();
+
+        for (int i = points.size()-1; i >= 0; --i) {
+            Point2d p = points.get(i);
+            spatialIndex.insert(p.envelope(), p);
+        }
+
+        int W = samples.width;
+        int H = samples.height;
+
+        double cellWidth  = boundingBox.getWidth()  / W;
+        double cellHeight = boundingBox.getHeight() / H;
+
+        Envelope     queryBuffer = new Envelope();
+        Point2d []   neighbors   = new Point2d[4];
+        Coordinate   center      = new Coordinate();
+        L1Comparator invL1       = new L1Comparator(center);
+
+        double [] raster = new double[W*H];
+        Arrays.fill(raster, Double.NaN);
+
+        int row = 0;
+        for (int j = 0; j < H; ++j, row += W) {
+            double y = j*cellHeight + 0.5d*cellHeight;
+            double x = 0.5d*cellWidth;
+            for (int i = 0; i < W; ++i, x += cellWidth) {
+                queryBuffer.init(
+                    x - dxMax, x + dxMax, 
+                    y - dyMax, y + dyMax);
+
+                List potential = spatialIndex.query(queryBuffer);
+
+                if (potential.isEmpty()) {
+                    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);
+                }
+            }
+        }
+
+        this.raster = raster;
+        this.width  = W;
+        this.height = H;
+
+        return true;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :