Mercurial > dive4elements > gnv-client
view gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java @ 592:3c939c95c477
Removed accidental checkin of Test-Shapefiles
gnv-artifacts/trunk@647 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Tim Englich <tim.englich@intevation.de> |
---|---|
date | Fri, 29 Jan 2010 13:28:28 +0000 |
parents | 44415ae01ddb |
children | 9a828e5a2390 |
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; /** * @author Sascha L. Teichmann (sascha.teichmann@intevation.de) */ public final class Interpolation2D { private static Logger log = Logger.getLogger(Interpolation2D.class); public static final int CULL_POINT_THRESHOLD = Integer.getInteger( "gnv.interpolation2d.cull.point.threshold", 1000); public static final double EPS = 1e-6d; public interface Consumer { void interpolated(Coordinate point, boolean success); } // interface Consumer private Interpolation2D() { } public static final Envelope relevantArea( List<? extends Coordinate> path, List<? extends Coordinate> points ) { return relevantArea(path, points, CULL_POINT_THRESHOLD); } public static final Envelope relevantArea( Envelope pathBBox, List<? extends Coordinate> points ) { return relevantArea(pathBBox, points, CULL_POINT_THRESHOLD); } public static final Envelope relevantArea( Envelope pathBBox, List<? extends Coordinate> points, int threshold ) { return points.size() < threshold ? null : relevantArea( pathBBox, pointsBoundingBox(points)); } 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)); } 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; } 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) }; } public static final double pos(double [] x) { return x[0] >= 0 ? x[0] : x[1]; } 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; } 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); } } 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: