--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java	Fri Sep 28 12:14:00 2012 +0200
@@ -0,0 +1,210 @@
+/*
+ * Copyright (c) 2010 by Intevation GmbH
+ *
+ * This program is free software under the LGPL (>=v2.1)
+ * Read the file LGPL.txt coming with the software for details
+ * or visit http://www.gnu.org/licenses/ if it does not exist.
+ */
+
+package de.intevation.gnv.math;
+
+import com.vividsolutions.jts.geom.Coordinate;
+import com.vividsolutions.jts.geom.Envelope;
+
+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;
+
+/**
+ * Does an area interpolation for "Horizontalschnitte".
+ *
+ * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
+ */
+public class AreaInterpolation
+implements   Serializable
+{
+    private static Logger log = Logger.getLogger(AreaInterpolation.class);
+
+    /**
+     * The generated raster.
+     */
+    protected double [] raster;
+
+    /**
+     * The width of the raster.
+     */
+    protected int width;
+    /**
+     * The height of the raster.
+     */
+    protected int height;
+
+    /**
+     * Default constructor.
+     */
+    public AreaInterpolation() {
+    }
+
+    /**
+     * Returns the width of the generated raster.
+     * @return the width of the raster.
+     */
+    public int getWidth() {
+        return width;
+    }
+
+    /**
+     * Returns the height of the raster.
+     * @return The height of the raster.
+     */
+    public int getHeight() {
+        return height;
+    }
+
+    /**
+     * The generated raster.
+     * @return the raster.
+     */
+    public double [] getRaster() {
+        return raster;
+    }
+
+    /**
+     * Epsilon for numerical stability.
+     */
+    public static final double EPS = 1e-6d;
+
+
+    /**
+     * Fills a raster by interpolating the input values.
+     * @param points The attributed input values.
+     * @param boundingBox The world area where to interpolate.
+     * @param samples The width and height of the raster.
+     * @param depth The callback to query the depth at a given point.
+     * @param extrapolationRounds Number of extrapolation point layers
+     * to generate before doing the interpolation.
+     * @return true if the interpolation succeeds else false.
+     */
+    public boolean interpolate(
+        List<? extends Point2d> points,
+        Envelope                boundingBox,
+        Dimension               samples,
+        XYDepth                 depth,
+        int                     extrapolationRounds
+    ) {
+        boolean debug = log.isDebugEnabled();
+
+        if (points == null || points.isEmpty()) {
+            log.warn("no points to interpolate");
+            return false;
+        }
+
+        List<GridCell> cells = GridCell.pointsToGridCells(
+            points,
+            Interpolation2D.relevantArea(
+                boundingBox,
+                points),
+            extrapolationRounds);
+
+        if (cells.isEmpty()) {
+            log.warn("no cells to interpolate");
+            return false;
+        }
+
+        int W = samples.width;
+        int H = samples.height;
+
+        double cellWidth  = boundingBox.getWidth()  / W;
+        double cellHeight = boundingBox.getHeight() / H;
+
+        STRtree spatialIndex = new STRtree();
+
+        for (GridCell cell: cells) {
+            spatialIndex.insert(cell.getEnvelope(), cell);
+        }
+
+        if (debug) {
+            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);
+        }
+
+        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);
+
+        double minX = boundingBox.getMinX();
+        double minY = boundingBox.getMinY();
+
+        long startTime = System.currentTimeMillis();
+
+        int pos = 0;
+        for (int j = 0; j < H; ++j) {
+
+            double y = j*cellHeight + 0.5d*cellHeight + minY;
+            double x = 0.5d*cellWidth + minX;
+
+            for (int end = pos + W; pos < end; ++pos, x += cellWidth) {
+
+                queryBuffer.init(x - EPS, x + EPS, y - EPS, y + EPS);
+                center.x = x; center.y = y;
+                finder.prepare(center);
+                spatialIndex.query(queryBuffer, finder);
+
+                GridCell found = finder.found;
+
+                if (found == null || depth.depth(center) > 0d) {
+                    continue;
+                }
+
+                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[pos] = Interpolation2D.interpolate(
+                    y1, z1,
+                    y2, z2,
+                    center.y);
+            }
+        }
+
+        long stopTime = System.currentTimeMillis();
+
+        if (debug) {
+            log.debug("interpolation took: " +
+                (stopTime - startTime)/1000f + " secs");
+        }
+
+        this.raster = raster;
+        this.width  = W;
+        this.height = H;
+
+        return true;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :