Mercurial > dive4elements > gnv-client
view gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java @ 1037:da4306683186
Added queries for TimeSeriesPoints in vectorvalues for TimeSeries
gnv-artifacts/trunk@1109 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Tim Englich <tim.englich@intevation.de> |
---|---|
date | Tue, 18 May 2010 14:27:46 +0000 |
parents | 05bf8534a35a |
children | f953c9a559d8 |
line wrap: on
line source
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 :