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 :