ingo@1115: /*
ingo@1115:  * Copyright (c) 2010 by Intevation GmbH
ingo@1115:  *
ingo@1115:  * This program is free software under the LGPL (>=v2.1)
ingo@1115:  * Read the file LGPL.txt coming with the software for details
ingo@1115:  * or visit http://www.gnu.org/licenses/ if it does not exist.
ingo@1115:  */
ingo@1115: 
sascha@361: package de.intevation.gnv.math;
sascha@361: 
sascha@361: import com.vividsolutions.jts.geom.Coordinate;
sascha@361: import com.vividsolutions.jts.geom.Envelope;
sascha@361: 
sascha@519: import com.vividsolutions.jts.index.strtree.STRtree;
sascha@501: 
sascha@501: import java.util.List;
sascha@501: 
ingo@365: import org.apache.log4j.Logger;
ingo@365: 
sascha@361: /**
sascha@807:  * Interpolates along a given line string. This used to generate
sascha@807:  * "horizontale Schnittprofile".
sascha@807:  *
sascha@807:  * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
sascha@361:  */
sascha@361: public final class Interpolation2D
sascha@361: {
ingo@365:     private static Logger log = Logger.getLogger(Interpolation2D.class);
ingo@365: 
sascha@807:     /**
sascha@807:      * If the number of input points are over this threshold
sascha@807:      * some culling strategies are applied to reduce the amount of points.
sascha@807:      */
sascha@519:     public static final int CULL_POINT_THRESHOLD = Integer.getInteger(
sascha@519:         "gnv.interpolation2d.cull.point.threshold", 1000);
sascha@519: 
sascha@807:     /**
sascha@807:      * Numerical stability epsilon.
sascha@807:      */
sascha@519:     public static final double EPS = 1e-6d;
sascha@519: 
sascha@807:     /**
sascha@807:      * Callback to send the interpolated values back to the
sascha@807:      * caller without building large temporary strcutures.
sascha@807:      */
sascha@361:     public interface Consumer {
sascha@807:         /**
sascha@807:          * Sends an interpolated point back to the called.
sascha@807:          * The interpolated parameter is stored in the z component.
sascha@807:          * @param point The interpolated point.
sascha@807:          * @param success true if interpolation at the point
sascha@807:          * succeed, else false.
sascha@807:          */
ingo@416:         void interpolated(Coordinate point, boolean success);
sascha@361:     } // interface Consumer
sascha@361: 
sascha@361:     private Interpolation2D() {
sascha@361:     }
sascha@361: 
sascha@807:     /**
sascha@807:      * Calculates the relevant area for a given line string.
sascha@807:      * @param path The line string.
sascha@807:      * @param points The sample points.
sascha@807:      * @return The relevant area.
sascha@807:      */
sascha@528:     public static final Envelope relevantArea(
sascha@519:         List<? extends Coordinate> path,
sascha@528:         List<? extends Coordinate> points
sascha@528:     ) {
sascha@528:         return relevantArea(path, points, CULL_POINT_THRESHOLD);
sascha@528:     }
sascha@528: 
sascha@807:     /**
sascha@807:      * Calculates the relevant area for a given bounding box.
sascha@807:      * @param pathBBox The bounding box.
sascha@807:      * @param points The sample points.
sascha@807:      * @return The relevant area.
sascha@807:      */
sascha@528:     public static final Envelope relevantArea(
sascha@528:         Envelope                   pathBBox,
sascha@528:         List<? extends Coordinate> points
sascha@528:     ) {
sascha@528:         return relevantArea(pathBBox, points, CULL_POINT_THRESHOLD);
sascha@528:     }
sascha@528: 
sascha@807:     /**
sascha@807:      * Calculates the relevant area for a given bounding box.
sascha@807:      * @param pathBBox The bounding box.
sascha@807:      * @param points The sample points.
sascha@807:      * @param threshold If the number of sample points
sascha@807:      * are below this threshold no relevant area is calculated.
sascha@807:      * @return The relevant area.
sascha@807:      */
sascha@528:     public static final Envelope relevantArea(
sascha@528:         Envelope                   pathBBox,
sascha@528:         List<? extends Coordinate> points,
sascha@528:         int                        threshold
sascha@528:     ) {
sascha@528:         return points.size() < threshold
sascha@528:             ? null
sascha@528:             : relevantArea(
sascha@528:                 pathBBox,
sascha@528:                 pointsBoundingBox(points));
sascha@528:     }
sascha@528: 
sascha@807:     /**
sascha@807:      * Calculates the relevant area for a given line string.
sascha@807:      * @param path The line string.
sascha@807:      * @param points The sample points.
sascha@807:      * @param threshold If the number of sample points
sascha@807:      * are below this threshold no relevant area is calculated.
sascha@807:      * @return The relevant area.
sascha@807:      */
sascha@528:     public static final Envelope relevantArea(
sascha@528:         List<? extends Coordinate> path,
sascha@528:         List<? extends Coordinate> points,
sascha@528:         int                        threshold
sascha@528:     ) {
sascha@528:         return points.size() < threshold || path.isEmpty()
sascha@528:             ? null
sascha@528:             : relevantArea(
sascha@528:                 pointsBoundingBox(path),
sascha@528:                 pointsBoundingBox(points));
sascha@528:     }
sascha@528: 
sascha@807:     /**
sascha@807:      * Heuristic function to define a 'relevant area'.
sascha@807:      * @param pathBBox The bounding box of the line line string.
sascha@807:      * @param pointsBBox The bounding box of the sample points.
sascha@807:      * @return The relevant area.
sascha@807:      */
sascha@528:     public static final Envelope relevantArea(
sascha@528:         Envelope pathBBox,
sascha@528:         Envelope pointsBBox
sascha@528:     ) {
sascha@528:         double pathArea   = pathBBox.getWidth()*pathBBox.getHeight();
sascha@528:         double pointsArea = pointsBBox.getWidth()*pointsBBox.getHeight();
sascha@528: 
sascha@528:         if (pathArea > 0.8d*pointsArea) { return null; }
sascha@528: 
sascha@528:         double nArea = 1.44d * pathArea;
sascha@528:         if (nArea < 0.1d*pointsArea) nArea = 0.1d*pointsArea;
sascha@528:         double w = pathBBox.getWidth();
sascha@528:         double h = pathBBox.getHeight();
sascha@528:         double [] d = solveQuadratic(1d, w+h, pathArea - nArea);
sascha@528: 
sascha@528:         if (d == null) { return null; }
sascha@528: 
sascha@528:         double extra = pos(d);
sascha@528: 
sascha@528:         pathBBox.expandBy(extra);
sascha@528: 
sascha@528:         return pathBBox;
sascha@528:     }
sascha@528: 
sascha@807:     /**
sascha@807:      * Solves quadratic equation a*x*x + b*x + c = 0.
sascha@807:      * @param a a-coefficent.
sascha@807:      * @param b b-coefficent
sascha@807:      * @param c c-coefficent.
sascha@807:      * @return the solution of the equation, or null if not solvable.
sascha@807:      */
sascha@528:     public static final double [] solveQuadratic(
sascha@528:         double a, double b, double c
sascha@528:     ) {
sascha@528:         double d = b*b - 4d*a*c;
sascha@528:         if (d < 0d) { return null; }
sascha@528: 
sascha@528:         d = Math.sqrt(d);
sascha@528:         a = 1d/(2d*a);
sascha@528:         b = -b;
sascha@528: 
sascha@528:         return new double [] { a*(b + d), a*(b - d) };
sascha@528:     }
sascha@528: 
sascha@807:     /**
sascha@807:      * Return the element of a two element array which
sascha@807:      * is greater or equal zero.
sascha@807:      * @param x The two values.
sascha@807:      * @return The value which is greater or equal zero.
sascha@807:      */
sascha@528:     public static final double pos(double [] x) {
sascha@528:         return x[0] >= 0 ? x[0] : x[1];
sascha@528:     }
sascha@528: 
sascha@528: 
sascha@807:     /**
sascha@807:      * Calculates the bounding box of a given line string.
sascha@807:      * @param path The line string.
sascha@807:      * @return The bounding box.
sascha@807:      */
sascha@528:     public static Envelope pointsBoundingBox(
sascha@528:         List<? extends Coordinate> path
sascha@501:     ) {
sascha@519:         int N = path.size();
sascha@519:         Envelope area = new Envelope(path.get(N-1));
sascha@501: 
sascha@519:         for (int i = N-2; i >= 0; --i) {
sascha@519:             area.expandToInclude(path.get(i));
sascha@501:         }
sascha@501: 
sascha@519:         return area;
sascha@431:     }
sascha@431: 
sascha@807:     /**
sascha@807:      * Interpolates linearly a number of coordinates and parameter values along
sascha@807:      * a given line string path. The results are issued to a consumer.
sascha@807:      * @param path The line string path.
sascha@807:      * @param points The sample points.
sascha@807:      * @param from Start point as a scalar value linear
sascha@807:      * referenced on the line string.
sascha@807:      * @param to End point of  as a scalar value linear
sascha@807:      * referenced on the line string.
sascha@807:      * @param steps Number of points to be interpolated.
sascha@807:      * @param metrics The used metric.
sascha@807:      * @param consumer The callback to retrieve the result points.
sascha@807:      */
sascha@361:     public static void interpolate(
sascha@361:         List<? extends Coordinate> path,
sascha@361:         List<? extends Point2d>    points,
sascha@361:         double                     from,
sascha@361:         double                     to,
sascha@361:         int                        steps,
sascha@361:         Metrics                    metrics,
sascha@361:         Consumer                   consumer
sascha@361:     ) {
sascha@434:         boolean debug = log.isDebugEnabled();
sascha@434: 
sascha@361:         int N = path.size();
sascha@361:         int M = points.size();
sascha@361: 
sascha@434:         if (debug) {
ingo@416:             log.debug("Size of path: " + N);
ingo@416:             log.debug("Size of points: " + M);
ingo@416:         }
ingo@365: 
sascha@361:         if (M < 1 || N < 2) { // nothing to do
sascha@361:             return;
sascha@361:         }
ingo@365: 
sascha@519:         List<GridCell> cells = GridCell.pointsToGridCells(
sascha@528:             points, relevantArea(path, points));
ingo@365: 
sascha@519:         if (cells.isEmpty()) {
sascha@519:             log.warn("no cells found");
sascha@519:             return;
ingo@416:         }
sascha@361: 
sascha@361:         // put into spatial index to speed up finding neighbors.
sascha@519:         STRtree spatialIndex = new STRtree();
sascha@361: 
sascha@519:         for (GridCell cell: cells) {
sascha@519:             spatialIndex.insert(cell.getEnvelope(), cell);
sascha@361:         }
sascha@361: 
sascha@361:         LinearToMap linearToMap = new LinearToMap(
sascha@361:             path, from, to, metrics);
sascha@361: 
sascha@361:         double dP = (to - from)/steps;
sascha@361: 
sascha@519:         Coordinate          center      = new Coordinate();
sascha@519:         Envelope            queryBuffer = new Envelope();
sascha@519:         GridCell.CellFinder finder      = new GridCell.CellFinder();
sascha@361: 
ingo@365:         int missedInterpolations = 0;
sascha@519:         int interpolations       = 0;
sascha@474: 
ingo@365:         for (double p = from; p <= to; p += dP) {
sascha@361:             if (!linearToMap.locate(p, center)) {
sascha@361:                 continue;
sascha@361:             }
sascha@361:             queryBuffer.init(
sascha@778:                 center.x - EPS, center.x + EPS,
sascha@519:                 center.y - EPS, center.y + EPS);
sascha@361: 
sascha@519:             finder.prepare(center);
sascha@519:             spatialIndex.query(queryBuffer, finder);
sascha@361: 
sascha@519:             GridCell found = finder.found;
sascha@519: 
sascha@519:             if (found == null) {
sascha@519:                 consumer.interpolated(center, false);
sascha@519:                 ++missedInterpolations;
sascha@519:                 continue;
sascha@361:             }
sascha@361: 
sascha@519:             Point2d n0 = found.p1;
sascha@519:             Point2d n1 = found.p2;
sascha@519:             Point2d n2 = found.p3;
sascha@519:             Point2d n3 = found.p4;
sascha@361: 
sascha@519:             double z1 = interpolate(
sascha@519:                 n0.x, n0.z,
sascha@519:                 n1.x, n1.z,
sascha@519:                 center.x);
sascha@519:             double z2 = interpolate(
sascha@519:                 n2.x, n2.z,
sascha@519:                 n3.x, n3.z,
sascha@519:                 center.x);
sascha@519:             double y1 = interpolate(
sascha@519:                 n0.x, n0.y,
sascha@519:                 n1.x, n1.y,
sascha@519:                 center.x);
sascha@519:             double y2 = interpolate(
sascha@519:                 n2.x, n2.y,
sascha@519:                 n3.x, n3.y,
sascha@519:                 center.x);
sascha@519:             center.z = interpolate(
sascha@519:                 y1, z1,
sascha@519:                 y2, z2,
sascha@519:                 center.y);
sascha@519:             consumer.interpolated(center, true);
sascha@519:             ++interpolations;
sascha@361:         }
ingo@365: 
sascha@434:         if (debug) {
sascha@778:             log.debug("interpolations: " +
ingo@416:                 interpolations + " / " + missedInterpolations);
ingo@416:         }
sascha@361:     }
sascha@361: 
sascha@807:     /**
sascha@807:      * Linear interpolate a value between (x1, y1) and (x2, y2) at
sascha@807:      * a given x-value.
sascha@807:      * @param x1 x component of first point.
sascha@807:      * @param y1 y component of first point.
sascha@807:      * @param x2 x component of second point.
sascha@807:      * @param y2 y component of second point.
sascha@807:      * @param x The x value.
sascha@807:      * @return The intepolated result.
sascha@807:      */
sascha@361:     public static final double interpolate(
sascha@361:         double x1, double y1,
sascha@361:         double x2, double y2,
sascha@361:         double x
sascha@361:     ) {
sascha@361:         if (x2 == x1) {
sascha@361:             return (y1 + y2)*0.5d;
sascha@361:         }
sascha@361:         double m = (y2-y1)/(x2-x1);
sascha@361:         double b = y1 - m*x1;
sascha@361:         return m*x + b;
sascha@361:     }
sascha@361: }
sascha@836: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :