Mercurial > dive4elements > gnv-client
diff gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java @ 361:aec85d00d82c
Added code to do 2D interpolations along a digitied track on the map.
Not connected to the transition model, yet.
gnv-artifacts/trunk@435 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 15 Dec 2009 22:25:53 +0000 |
parents | |
children | f66088a43ecc |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java Tue Dec 15 22:25:53 2009 +0000 @@ -0,0 +1,161 @@ +package de.intevation.gnv.math; + +import java.util.List; +import java.util.Collections; + +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.Envelope; + +import com.vividsolutions.jts.index.quadtree.Quadtree; + +/** + * @author Sascha L. Teichmann + */ +public final class Interpolation2D +{ + public interface Consumer { + void interpolated(Coordinate point); + } // interface Consumer + + private Interpolation2D() { + } + + public static void interpolate( + List<? extends Coordinate> path, + List<? extends Point2d> points, + double from, + double to, + int steps, + Metrics metrics, + Consumer consumer + ) { + int N = path.size(); + int M = points.size(); + + if (M < 1 || N < 2) { // nothing to do + return; + } + // figure out max delta(p[i].x, p[i-1].x) + Collections.sort(points, Point2d.X_COMPARATOR); + double dxMax = -Double.MAX_VALUE; + for (int i = 1; i < M; ++i) { + double dx = Math.abs(path.get(i).x - path.get(i-1).x); + if (dx > dxMax) { + dxMax = dx; + } + } + + dxMax = dxMax*0.5d + 1e-5d; + + // figure out max delta(p[i].y, p[i-1].y) + Collections.sort(path, Point2d.X_COMPARATOR); + double dyMax = -Double.MAX_VALUE; + for (int i = 1; i < M; ++i) { + double dy = Math.abs(path.get(i).y - path.get(i-1).y); + if (dy > dyMax) { + dyMax = dy; + } + } + + dyMax = dyMax*0.5d + 1e-5d; + + // put into spatial index to speed up finding neighbors. + Quadtree spatialIndex = new Quadtree(); + + for (int i = 0; i < M; ++i) { + Point2d p = points.get(i); + spatialIndex.insert(p.envelope(), p); + } + + LinearToMap linearToMap = new LinearToMap( + path, from, to, metrics); + + double dP = (to - from)/steps; + + Coordinate center = new Coordinate(); + + Envelope queryBuffer = new Envelope(); + + Point2d [] neighbors = new Point2d[4]; + + for (double p = to; p <= from; p += dP) { + if (!linearToMap.locate(p, center)) { + continue; + } + queryBuffer.init( + center.x - dxMax, center.x + dxMax, + center.y - dyMax, center.y + dyMax); + + List potential = spatialIndex.query(queryBuffer); + + L1Comparator invL1 = new L1Comparator(center); + Collections.sort(potential, invL1); + + neighbors[0] = neighbors[1] = + neighbors[2] = neighbors[3] = null; + + /* bit code of neighbors + 0---1 + | x | + 2---3 + */ + + int mask = 0; + // reversed order is essential here! + for (int i = potential.size()-1; i >= 0; --i) { + Point2d n = (Point2d)potential.get(i); + int code = n.x > center.x ? 1 : 0; + if (n.y > center.y) code |= 2; + neighbors[code] = n; + mask |= 1 << code; + } + + int numNeighbors = Integer.bitCount(mask); + + // only interpolate if we have all four neighbors + // and we do not have any gaps. + if (numNeighbors == 4 + && !neighbors[0].hasIGap(neighbors[0]) + && !neighbors[1].hasJGap(neighbors[3]) + && !neighbors[3].hasIGap(neighbors[2]) + && !neighbors[2].hasJGap(neighbors[0]) + ) { + double z1 = interpolate( + neighbors[0].x, neighbors[0].z, + neighbors[1].x, neighbors[1].z, + center.x); + double z2 = interpolate( + neighbors[2].x, neighbors[2].z, + neighbors[3].x, neighbors[3].z, + center.x); + double y1 = interpolate( + neighbors[0].x, neighbors[0].y, + neighbors[1].x, neighbors[1].y, + center.x); + double y2 = interpolate( + neighbors[2].x, neighbors[2].y, + neighbors[3].x, neighbors[3].y, + center.x); + center.z = interpolate( + y1, z1, + y2, z2, + center.y); + consumer.interpolated(center); + } + } + } + + 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: