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<? 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@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@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@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@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<? 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@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@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: