--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java	Fri Sep 28 12:13:56 2012 +0200
@@ -0,0 +1,343 @@
+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;
+
+/**
+ * Interpolates parameter values along a given line string from surface
+ * to the sea ground to generate "Profilschnitte".
+ *
+ * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
+ */
+public class Interpolation3D
+implements   Serializable
+{
+    private static Logger log = Logger.getLogger(Interpolation3D.class);
+
+    /**
+     * The default width of the interpolation: {@value}
+     */
+    public static final int DEFAULT_WIDTH  = 1024;
+
+    /**
+     * The default height of the interpolation: {@value}
+     */
+    public static final int DEFAULT_HEIGHT =  768;
+
+    /**
+     * Epsilon for numerical stability.
+     */
+    public static final double EPS = 1e-6d;
+
+    /**
+     * The width of the interpolation.
+     */
+    protected int width;
+
+    /**
+     * The height of the interpolation.
+     */
+    protected int height;
+
+    /**
+     * The cell width of the interpolation in world units.
+     */
+    protected double cellWidth;
+
+    /**
+     * The cell height of the interpolation in world units.
+     */
+    protected double cellHeight;
+
+    /**
+     * The coordinates of the interpolation.
+     */
+    protected Coordinate[] coordinates;
+
+    /**
+     * The generated raster.
+     */
+    protected double [] raster;
+
+    /**
+     * The sea ground depth along the line string.
+     */
+    protected double [] depths;
+
+    /**
+     * Default constructor.
+     */
+    public Interpolation3D() {
+        this(DEFAULT_WIDTH, DEFAULT_HEIGHT);
+    }
+
+    /**
+     * Constructor to create a Interpolation3D with a given size.
+     * @param size The size of the interpolation.
+     */
+    public Interpolation3D(Dimension size) {
+        this(size.width, size.height);
+    }
+
+    /**
+     * Constructor to create a Interpolation3D with a given size.
+     * @param width The width of the interpolation.
+     * @param height the height of the interpolation.
+     */
+    public Interpolation3D(int width, int height) {
+        this.width  = width;
+        this.height = height;
+    }
+
+    /**
+     * Returns the raster height of the interpolation.
+     * @return The raster height of the interpolation.
+     */
+    public int getHeight() {
+        return height;
+    }
+
+    /**
+     * Returns the raster width of the interpolation.
+     * @return The raster width of the interpolation.
+     */
+    public int getWidth() {
+        return width;
+    }
+
+    /**
+     * Returns the cell width of the interpolation in world units.
+     * @return The cell width of the interpolation in world units.
+     */
+    public double getCellWidth() {
+        return cellWidth;
+    }
+
+    /**
+     * Returns the cell height of the interpolation in world units.
+     * @return The cell height of the interpolation in world units.
+     */
+    public double getCellHeight() {
+        return cellHeight;
+    }
+
+    /**
+     * Returns the coordinates used for the interpolation.
+     * @return the coordinates used for the interpolation.
+     */
+    public Coordinate[] getCoordinates() {
+        return coordinates;
+    }
+
+    /**
+     * Returns the generated raster.
+     * @return The generated raster.
+     */
+    public double [] getRaster() {
+        return raster;
+    }
+
+    /**
+     * Returns the depths along the line string path.
+     * @return The depth along the line string path.
+     */
+    public double [] getDepths() {
+        return depths;
+    }
+
+    /**
+     * Returns the deepest depth along the line string path.
+     * @return The deepest depth along the line string path.
+     */
+    public double getMaxDepth() {
+        double maxDepth = Double.MAX_VALUE;
+        for (int i = depths!=null?depths.length-1:0; i >= 0; --i) {
+            double d = depths[i];
+            if (!Double.isNaN(d) && d < maxDepth) {
+                maxDepth = d;
+            }
+        }
+        return maxDepth;
+    }
+
+    /**
+     * Interpolates parameters along a given line string path from the surface
+     * to the sea ground.
+     * @param path The line string path.
+     * @param points The sample points.
+     * @param from Start point in scalar terms of the line string.
+     * @param to End point in scalar terms of the line string.
+     * @param metrics The used metric.
+     * @param xyDepth The callback to query the depth at a given point.
+     * @return true if the interpolation succeeds, else false.
+     */
+    public boolean interpolate(
+        List<? extends Coordinate> path,
+        List<? extends XYColumn>   points,
+        double                     from,
+        double                     to,
+        Metrics                    metrics,
+        XYDepth                    xyDepth
+    ) {
+        boolean debug = log.isDebugEnabled();
+
+        int N = path.size();
+        int M = points.size();
+
+        if (debug) {
+            log.debug("Size of path: " + N);
+            log.debug("Size of points: " + M);
+        }
+
+        if (M < 1 || N < 2) { // nothing to do
+            return false;
+        }
+
+        List<GridCell> cells = GridCell.pointsToGridCells(
+            points, Interpolation2D.relevantArea(path, points));
+
+        if (cells.isEmpty()) {
+            log.warn("no cells found");
+            return false;
+        }
+
+        // put into spatial index to speed up finding neighbors.
+        STRtree spatialIndex = new STRtree();
+
+        for (GridCell cell: cells) {
+            spatialIndex.insert(cell.getEnvelope(), cell);
+        }
+
+        LinearToMap linearToMap = new LinearToMap(
+            path, from, to, metrics);
+
+        double [] depths = new double[width];
+        Arrays.fill(depths, Double.NaN);
+
+        Coordinate[] coordinates = new Coordinate[width];
+
+        double cellWidth = (to - from)/width;
+
+        double maxDepth = Double.MAX_VALUE;
+
+        int i = 0;
+        Coordinate center = new Coordinate();
+        for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
+            if (linearToMap.locate(p, center)) {
+                coordinates[i] = (Coordinate) center.clone();
+                double depth = xyDepth.depth(center);
+                depths[i] = depth;
+                if (depth < maxDepth) {
+                    maxDepth = depth;
+                }
+            }
+        }
+
+        if (maxDepth == Double.MAX_VALUE) {
+            log.warn("no depth found -> no interpolation");
+            return false;
+        }
+
+        double cellHeight = Math.abs(maxDepth)/height;
+
+        if (debug) {
+            log.debug("max depth found: " + maxDepth);
+            log.debug("cell size: " + cellWidth + " x " + cellHeight);
+        }
+
+        double [] raster = new double[width*height];
+        Arrays.fill(raster, Double.NaN);
+
+        Envelope            queryBuffer  = new Envelope();
+        GridCell.CellFinder finder       = new GridCell.CellFinder();
+
+        i = 0;
+        for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
+            double depth = depths[i];
+            if (Double.isNaN(depth) || depth >= 0d) {
+                continue;
+            }
+            linearToMap.locate(p, center);
+
+            queryBuffer.init(
+                center.x - EPS, center.x + EPS,
+                center.y - EPS, center.y + EPS);
+
+            finder.prepare(center);
+            spatialIndex.query(queryBuffer, finder);
+
+            GridCell found = finder.found;
+
+            if (found == null) {
+                continue;
+            }
+
+            XYColumn n0 = (XYColumn)found.p1;
+            XYColumn n1 = (XYColumn)found.p2;
+            XYColumn n2 = (XYColumn)found.p3;
+            XYColumn n3 = (XYColumn)found.p4;
+
+            if (n0.prepare(xyDepth)
+            &&  n1.prepare(xyDepth)
+            &&  n2.prepare(xyDepth)
+            &&  n3.prepare(xyDepth)
+            ) {
+                double y1 = Interpolation2D.interpolate(
+                    n0.x, n0.y,
+                    n1.x, n1.y,
+                    center.x);
+                double y2 = Interpolation2D.interpolate(
+                    n2.x, n2.y,
+                    n3.x, n3.y,
+                    center.x);
+                double z = -cellHeight*0.5;
+                for (int j = i;
+                    j < raster.length && z >= depth;
+                    z -= cellHeight, j += width) {
+
+                    double v0 = n0.value(z);
+                    double v1 = n1.value(z);
+                    double v2 = n2.value(z);
+                    double v3 = n3.value(z);
+
+                    double z1 = Interpolation2D.interpolate(
+                        n0.x, v0,
+                        n1.x, v1,
+                        center.x);
+                    double z2 = Interpolation2D.interpolate(
+                        n2.x, v2,
+                        n3.x, v3,
+                        center.x);
+                    double value = Interpolation2D.interpolate(
+                        y1, z1,
+                        y2, z2,
+                        center.y);
+                    raster[j] = value;
+                }
+                // XXX: Adjusted depth to create no gap
+                // between last value and ground.
+                depths[i] = z+0.5d*cellHeight;
+            } // down the x/y column
+        } // all along the track
+
+        this.coordinates = coordinates;
+        this.depths      = depths;
+        this.raster      = raster;
+        this.cellWidth   = cellWidth;
+        this.cellHeight  = cellHeight;
+
+        return true;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :