--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Fri Sep 28 12:13:56 2012 +0200
@@ -0,0 +1,337 @@
+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.util.List;
+
+import org.apache.log4j.Logger;
+
+/**
+ * Interpolates along a given line string. This used to generate
+ * "horizontale Schnittprofile".
+ *
+ * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
+ */
+public final class Interpolation2D
+{
+    private static Logger log = Logger.getLogger(Interpolation2D.class);
+
+    /**
+     * If the number of input points are over this threshold
+     * some culling strategies are applied to reduce the amount of points.
+     */
+    public static final int CULL_POINT_THRESHOLD = Integer.getInteger(
+        "gnv.interpolation2d.cull.point.threshold", 1000);
+
+    /**
+     * Numerical stability epsilon.
+     */
+    public static final double EPS = 1e-6d;
+
+    /**
+     * Callback to send the interpolated values back to the
+     * caller without building large temporary strcutures.
+     */
+    public interface Consumer {
+        /**
+         * Sends an interpolated point back to the called.
+         * The interpolated parameter is stored in the z component.
+         * @param point The interpolated point.
+         * @param success true if interpolation at the point
+         * succeed, else false.
+         */
+        void interpolated(Coordinate point, boolean success);
+    } // interface Consumer
+
+    private Interpolation2D() {
+    }
+
+    /**
+     * Calculates the relevant area for a given line string.
+     * @param path The line string.
+     * @param points The sample points.
+     * @return The relevant area.
+     */
+    public static final Envelope relevantArea(
+        List<? extends Coordinate> path,
+        List<? extends Coordinate> points
+    ) {
+        return relevantArea(path, points, CULL_POINT_THRESHOLD);
+    }
+
+    /**
+     * Calculates the relevant area for a given bounding box.
+     * @param pathBBox The bounding box.
+     * @param points The sample points.
+     * @return The relevant area.
+     */
+    public static final Envelope relevantArea(
+        Envelope                   pathBBox,
+        List<? extends Coordinate> points
+    ) {
+        return relevantArea(pathBBox, points, CULL_POINT_THRESHOLD);
+    }
+
+    /**
+     * Calculates the relevant area for a given bounding box.
+     * @param pathBBox The bounding box.
+     * @param points The sample points.
+     * @param threshold If the number of sample points
+     * are below this threshold no relevant area is calculated.
+     * @return The relevant area.
+     */
+    public static final Envelope relevantArea(
+        Envelope                   pathBBox,
+        List<? extends Coordinate> points,
+        int                        threshold
+    ) {
+        return points.size() < threshold
+            ? null
+            : relevantArea(
+                pathBBox,
+                pointsBoundingBox(points));
+    }
+
+    /**
+     * Calculates the relevant area for a given line string.
+     * @param path The line string.
+     * @param points The sample points.
+     * @param threshold If the number of sample points
+     * are below this threshold no relevant area is calculated.
+     * @return The relevant area.
+     */
+    public static final Envelope relevantArea(
+        List<? extends Coordinate> path,
+        List<? extends Coordinate> points,
+        int                        threshold
+    ) {
+        return points.size() < threshold || path.isEmpty()
+            ? null
+            : relevantArea(
+                pointsBoundingBox(path),
+                pointsBoundingBox(points));
+    }
+
+    /**
+     * Heuristic function to define a 'relevant area'.
+     * @param pathBBox The bounding box of the line line string.
+     * @param pointsBBox The bounding box of the sample points.
+     * @return The relevant area.
+     */
+    public static final Envelope relevantArea(
+        Envelope pathBBox,
+        Envelope pointsBBox
+    ) {
+        double pathArea   = pathBBox.getWidth()*pathBBox.getHeight();
+        double pointsArea = pointsBBox.getWidth()*pointsBBox.getHeight();
+
+        if (pathArea > 0.8d*pointsArea) { return null; }
+
+        double nArea = 1.44d * pathArea;
+        if (nArea < 0.1d*pointsArea) nArea = 0.1d*pointsArea;
+        double w = pathBBox.getWidth();
+        double h = pathBBox.getHeight();
+        double [] d = solveQuadratic(1d, w+h, pathArea - nArea);
+
+        if (d == null) { return null; }
+
+        double extra = pos(d);
+
+        pathBBox.expandBy(extra);
+
+        return pathBBox;
+    }
+
+    /**
+     * Solves quadratic equation a*x*x + b*x + c = 0.
+     * @param a a-coefficent.
+     * @param b b-coefficent
+     * @param c c-coefficent.
+     * @return the solution of the equation, or null if not solvable.
+     */
+    public static final double [] solveQuadratic(
+        double a, double b, double c
+    ) {
+        double d = b*b - 4d*a*c;
+        if (d < 0d) { return null; }
+
+        d = Math.sqrt(d);
+        a = 1d/(2d*a);
+        b = -b;
+
+        return new double [] { a*(b + d), a*(b - d) };
+    }
+
+    /**
+     * Return the element of a two element array which
+     * is greater or equal zero.
+     * @param x The two values.
+     * @return The value which is greater or equal zero.
+     */
+    public static final double pos(double [] x) {
+        return x[0] >= 0 ? x[0] : x[1];
+    }
+
+
+    /**
+     * Calculates the bounding box of a given line string.
+     * @param path The line string.
+     * @return The bounding box.
+     */
+    public static Envelope pointsBoundingBox(
+        List<? extends Coordinate> path
+    ) {
+        int N = path.size();
+        Envelope area = new Envelope(path.get(N-1));
+
+        for (int i = N-2; i >= 0; --i) {
+            area.expandToInclude(path.get(i));
+        }
+
+        return area;
+    }
+
+    /**
+     * Interpolates linearly a number of coordinates and parameter values along
+     * a given line string path. The results are issued to a consumer.
+     * @param path The line string path.
+     * @param points The sample points.
+     * @param from Start point as a scalar value linear
+     * referenced on the line string.
+     * @param to End point of  as a scalar value linear
+     * referenced on the line string.
+     * @param steps Number of points to be interpolated.
+     * @param metrics The used metric.
+     * @param consumer The callback to retrieve the result points.
+     */
+    public static void interpolate(
+        List<? extends Coordinate> path,
+        List<? extends Point2d>    points,
+        double                     from,
+        double                     to,
+        int                        steps,
+        Metrics                    metrics,
+        Consumer                   consumer
+    ) {
+        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;
+        }
+
+        List<GridCell> cells = GridCell.pointsToGridCells(
+            points, relevantArea(path, points));
+
+        if (cells.isEmpty()) {
+            log.warn("no cells found");
+            return;
+        }
+
+        // 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 dP = (to - from)/steps;
+
+        Coordinate          center      = new Coordinate();
+        Envelope            queryBuffer = new Envelope();
+        GridCell.CellFinder finder      = new GridCell.CellFinder();
+
+        int missedInterpolations = 0;
+        int interpolations       = 0;
+
+        for (double p = from; p <= to; p += dP) {
+            if (!linearToMap.locate(p, center)) {
+                continue;
+            }
+            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) {
+                consumer.interpolated(center, false);
+                ++missedInterpolations;
+                continue;
+            }
+
+            Point2d n0 = found.p1;
+            Point2d n1 = found.p2;
+            Point2d n2 = found.p3;
+            Point2d n3 = found.p4;
+
+            double z1 = interpolate(
+                n0.x, n0.z,
+                n1.x, n1.z,
+                center.x);
+            double z2 = interpolate(
+                n2.x, n2.z,
+                n3.x, n3.z,
+                center.x);
+            double y1 = interpolate(
+                n0.x, n0.y,
+                n1.x, n1.y,
+                center.x);
+            double y2 = interpolate(
+                n2.x, n2.y,
+                n3.x, n3.y,
+                center.x);
+            center.z = interpolate(
+                y1, z1,
+                y2, z2,
+                center.y);
+            consumer.interpolated(center, true);
+            ++interpolations;
+        }
+
+        if (debug) {
+            log.debug("interpolations: " +
+                interpolations + " / " + missedInterpolations);
+        }
+    }
+
+    /**
+     * Linear interpolate a value between (x1, y1) and (x2, y2) at
+     * a given x-value.
+     * @param x1 x component of first point.
+     * @param y1 y component of first point.
+     * @param x2 x component of second point.
+     * @param y2 y component of second point.
+     * @param x The x value.
+     * @return The intepolated result.
+     */
+    public static final double interpolate(
+        double x1, double y1,
+        double x2, double y2,
+        double x
+    ) {
+        if (x2 == x1) {
+            return (y1 + y2)*0.5d;
+        }
+        double m = (y2-y1)/(x2-x1);
+        double b = y1 - m*x1;
+        return m*x + b;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :