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@780:  *  @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
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@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@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@798: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :