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 Sascha L. Teichmann 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 path, sascha@528: List 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 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 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 path, sascha@528: List 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 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 path, sascha@361: List 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 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 :