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@501: * @author Sascha L. Teichmann (sascha.teichmann@intevation.de) sascha@361: */ sascha@361: public final class Interpolation2D sascha@361: { ingo@365: private static Logger log = Logger.getLogger(Interpolation2D.class); ingo@365: sascha@519: public static final int CULL_POINT_THRESHOLD = Integer.getInteger( sascha@519: "gnv.interpolation2d.cull.point.threshold", 1000); sascha@519: sascha@519: public static final double EPS = 1e-6d; sascha@519: sascha@361: public interface Consumer { ingo@416: void interpolated(Coordinate point, boolean success); sascha@361: } // interface Consumer sascha@361: sascha@361: private Interpolation2D() { sascha@361: } sascha@361: 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@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@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@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@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@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@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@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@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@519: 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) { ingo@416: log.debug("interpolations: " + ingo@416: interpolations + " / " + missedInterpolations); ingo@416: } sascha@361: } sascha@361: 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@361: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8: