sascha@361: package de.intevation.gnv.math;
sascha@361:
ingo@365: import java.util.ArrayList;
sascha@361: import java.util.List;
ingo@365: import java.util.HashMap;
sascha@361: import java.util.Collections;
sascha@361:
sascha@361: import com.vividsolutions.jts.geom.Coordinate;
sascha@361: import com.vividsolutions.jts.geom.Envelope;
sascha@361:
sascha@361: import com.vividsolutions.jts.index.quadtree.Quadtree;
sascha@361:
ingo@365: import org.apache.log4j.Logger;
ingo@365:
sascha@361: /**
sascha@361: * @author Sascha L. Teichmann
sascha@361: */
sascha@361: public final class Interpolation2D
sascha@361: {
ingo@365: private static Logger log = Logger.getLogger(Interpolation2D.class);
ingo@365:
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@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@361: int N = path.size();
sascha@361: int M = points.size();
sascha@361:
ingo@416: if (log.isDebugEnabled()) {
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:
ingo@365: HashMap<Integer, ArrayList<Point2d>> map = new HashMap<Integer, ArrayList<Point2d>>();
ingo@365:
ingo@365: for (int k = M-1; k >= 0; --k) {
ingo@365: Point2d p = points.get(k);
ingo@365:
ingo@365: ArrayList<Point2d> list = map.get(p.j);
ingo@365:
ingo@365: if (list == null) {
ingo@365: map.put(p.j, list = new ArrayList<Point2d>());
ingo@365: }
ingo@365: list.add(p);
ingo@365: }
ingo@365:
sascha@361: double dxMax = -Double.MAX_VALUE;
ingo@365:
ingo@365: for (ArrayList<Point2d> v: map.values()) {
ingo@365: Collections.sort(v, Point2d.X_COMPARATOR);
ingo@365: for (int i = 1, L = v.size(); i < L; ++i) {
ingo@365: double dx = Math.abs(v.get(i).x - v.get(i-1).x);
ingo@365: if (dx > dxMax) {
ingo@365: dxMax = dx;
ingo@365: }
sascha@361: }
sascha@361: }
sascha@361:
ingo@416: dxMax += 1e-5d;
sascha@361:
ingo@365: map.clear();
ingo@365:
ingo@365: for (int k = M-1; k >= 0; --k) {
ingo@365: Point2d p = points.get(k);
ingo@365:
ingo@365: ArrayList<Point2d> list = map.get(p.i);
ingo@365:
ingo@365: if (list == null) {
ingo@365: map.put(p.i, list = new ArrayList<Point2d>());
ingo@365: }
ingo@365: list.add(p);
ingo@365: }
ingo@365:
sascha@361: double dyMax = -Double.MAX_VALUE;
ingo@365:
ingo@365: for (ArrayList<Point2d> v: map.values()) {
ingo@365: Collections.sort(v, Point2d.Y_COMPARATOR);
ingo@365: for (int i = 1, L = v.size(); i < L; ++i) {
ingo@365: double dy = Math.abs(v.get(i).y - v.get(i-1).y);
ingo@365: if (dy > dyMax) {
ingo@365: dyMax = dy;
ingo@365: }
sascha@361: }
sascha@361: }
sascha@361:
ingo@416: dyMax += 1e-5d;
ingo@365:
ingo@365: map = null;
ingo@365:
ingo@416: if (log.isDebugEnabled()) {
ingo@416: log.debug("buffer size: " + dxMax + " / " + dyMax);
ingo@416: }
sascha@361:
sascha@361: // put into spatial index to speed up finding neighbors.
sascha@361: Quadtree spatialIndex = new Quadtree();
sascha@361:
sascha@361: for (int i = 0; i < M; ++i) {
sascha@361: Point2d p = points.get(i);
sascha@361: spatialIndex.insert(p.envelope(), p);
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@361: Coordinate center = new Coordinate();
sascha@361:
sascha@361: Envelope queryBuffer = new Envelope();
sascha@361:
sascha@361: Point2d [] neighbors = new Point2d[4];
sascha@361:
ingo@365: int missedInterpolations = 0;
ingo@365: int interpolations = 0;
ingo@365:
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@361: center.x - dxMax, center.x + dxMax,
sascha@361: center.y - dyMax, center.y + dyMax);
sascha@361:
sascha@361: List potential = spatialIndex.query(queryBuffer);
sascha@361:
sascha@361: L1Comparator invL1 = new L1Comparator(center);
sascha@361: Collections.sort(potential, invL1);
sascha@361:
sascha@361: neighbors[0] = neighbors[1] =
sascha@361: neighbors[2] = neighbors[3] = null;
sascha@361:
sascha@361: /* bit code of neighbors
sascha@361: 0---1
sascha@361: | x |
sascha@361: 2---3
sascha@361: */
sascha@361:
sascha@361: int mask = 0;
sascha@361: // reversed order is essential here!
sascha@361: for (int i = potential.size()-1; i >= 0; --i) {
sascha@361: Point2d n = (Point2d)potential.get(i);
sascha@361: int code = n.x > center.x ? 1 : 0;
sascha@361: if (n.y > center.y) code |= 2;
sascha@361: neighbors[code] = n;
sascha@361: mask |= 1 << code;
sascha@361: }
sascha@361:
sascha@361: int numNeighbors = Integer.bitCount(mask);
sascha@361:
sascha@361: // only interpolate if we have all four neighbors
sascha@361: // and we do not have any gaps.
sascha@361: if (numNeighbors == 4
sascha@366: && !neighbors[0].hasIGap(neighbors[1])
sascha@361: && !neighbors[1].hasJGap(neighbors[3])
sascha@361: && !neighbors[3].hasIGap(neighbors[2])
sascha@361: && !neighbors[2].hasJGap(neighbors[0])
sascha@361: ) {
sascha@361: double z1 = interpolate(
sascha@361: neighbors[0].x, neighbors[0].z,
sascha@361: neighbors[1].x, neighbors[1].z,
sascha@361: center.x);
sascha@361: double z2 = interpolate(
sascha@361: neighbors[2].x, neighbors[2].z,
sascha@361: neighbors[3].x, neighbors[3].z,
sascha@361: center.x);
sascha@361: double y1 = interpolate(
sascha@361: neighbors[0].x, neighbors[0].y,
sascha@361: neighbors[1].x, neighbors[1].y,
sascha@361: center.x);
sascha@361: double y2 = interpolate(
sascha@361: neighbors[2].x, neighbors[2].y,
sascha@361: neighbors[3].x, neighbors[3].y,
sascha@361: center.x);
sascha@361: center.z = interpolate(
sascha@361: y1, z1,
sascha@361: y2, z2,
sascha@361: center.y);
ingo@416: consumer.interpolated(center, true);
ingo@365: ++interpolations;
ingo@365: }
ingo@365: else {
ingo@416: consumer.interpolated(center, false);
ingo@365: ++missedInterpolations;
sascha@361: }
sascha@361: }
ingo@365:
ingo@416: if (log.isDebugEnabled()) {
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: