# HG changeset patch # User Sascha L. Teichmann # Date 1260915953 0 # Node ID aec85d00d82caf0a73296c83ea614d2f1330960e # Parent ee760729f6b8a2771b65618d2a71e59947086d15 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 diff -r ee760729f6b8 -r aec85d00d82c gnv-artifacts/ChangeLog --- a/gnv-artifacts/ChangeLog Tue Dec 15 15:58:58 2009 +0000 +++ b/gnv-artifacts/ChangeLog Tue Dec 15 22:25:53 2009 +0000 @@ -1,3 +1,41 @@ +2009-12-15 Sascha L. Teichmann + + * src/main/java/de/intevation/gnv/math/LinearToMap.java: + Uses JTS Coordinate as geometry model now. + + * src/main/java/de/intevation/gnv/math/Metrics.java, + src/main/java/de/intevation/gnv/math/Interpolator.java: New. + Moved from inner class of LinearToMap to top level class + to be more reusable. Uses JTS Coordinate as geometry model now. + + * src/main/java/de/intevation/gnv/math/Point2d.java: New. + Extends JTS Coordinate to have an additional (i, j) + to model the topological neighborhood withi the mesh, too. + + * src/main/java/de/intevation/gnv/math/Interpolation2D.java: New. + Has a method interpolate() which takes a path line string in form + of a list of JTS Coordinates, a list of grid points (Point2d + to carry the topology, too), a linear range in diagram coordinate + space, a metric to cope with the projection. It reports + interpolated points to an implementor of the new inner interface + Consumer as a JTS Coordinate. (x, y) of this coordinate is the + postion on the map, the z value is the interpolated attribute. + + To speed up the search for the neighbors the input points are + sorted into a quadtree and are queried first level with a buffer of + size (max(abs(p[i].x - p[i+1].x)), max(abs(p[i].y - p[i+1].y))) + around the point to be interpolated. The second level filter + is performed by an inverse L1-ordering with region coding, so + that only the nearest four neighbors are taken into acount. + Only if all four neighbors are present and no + i- or j-gaps exist the interpolation is performed. TODO: Create + a better extrapolation strategy in these cases were these conditions + are not fulfilled. + + * src/main/java/de/intevation/gnv/state/profile/horizontal/HorizontalProfileMeshCrossOutputState.java: + Added a process() method to perform the interpolation. It does + nothing by now. TODO: bring it to life. + 2009-12-15 Sascha L. Teichmann * src/main/java/de/intevation/gnv/math/LinearToMap.java: Map linear diff -r ee760729f6b8 -r aec85d00d82c gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java --- /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 path, + List 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: diff -r ee760729f6b8 -r aec85d00d82c gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolator.java --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolator.java Tue Dec 15 22:25:53 2009 +0000 @@ -0,0 +1,12 @@ +package de.intevation.gnv.math; + +import com.vividsolutions.jts.geom.Coordinate; + +/** + * @author Sascha L. Teichmann + */ +public interface Interpolator +{ + void interpolate(double t, Coordinate v); +} +// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8: diff -r ee760729f6b8 -r aec85d00d82c gnv-artifacts/src/main/java/de/intevation/gnv/math/L1Comparator.java --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/L1Comparator.java Tue Dec 15 22:25:53 2009 +0000 @@ -0,0 +1,30 @@ +package de.intevation.gnv.math; + +import java.util.Comparator; + +import com.vividsolutions.jts.geom.Coordinate; + +/** + * @author Sascha L. Teichmann + */ +public class L1Comparator +implements Comparator +{ + private Coordinate ref; + + public L1Comparator(Coordinate ref) { + this.ref = ref; + } + + public int compare(Object a, Object b) { + Coordinate pa = (Coordinate)a; + Coordinate pb = (Coordinate)b; + double da = Point2d.L1(ref, pa); + double db = Point2d.L1(ref, pb); + if (da < db) return -1; + if (da > db) return +1; + return 0; + } +} +// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8: + diff -r ee760729f6b8 -r aec85d00d82c gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearMetrics.java --- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearMetrics.java Tue Dec 15 15:58:58 2009 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearMetrics.java Tue Dec 15 22:25:53 2009 +0000 @@ -1,9 +1,6 @@ package de.intevation.gnv.math; -import java.awt.geom.Point2D; - -import de.intevation.gnv.math.LinearToMap.Metrics; -import de.intevation.gnv.math.LinearToMap.Interpolator; +import com.vividsolutions.jts.geom.Coordinate; /** * @author Sascha L. Teichmann @@ -22,7 +19,7 @@ private double my; private double by; - public LinearInterpolator(Point2D p1, Point2D p2) { + public LinearInterpolator(Coordinate p1, Coordinate p2) { /* I) p1.x = 0*m + bx @@ -35,28 +32,27 @@ mx = p2.x - p1.x */ - bx = p1.getX(); - mx = p2.getX() - bx; + bx = p1.x; + mx = p2.x - bx; - by = p1.getY(); - my = p2.getY() - by; + by = p1.y; + my = p2.y - by; } - public void interpolate(double t, Point2D v) { - double x = t*mx + bx; - double y = t*my + by; - v.setLocation(x, y); + public void interpolate(double t, Coordinate v) { + v.x = t*mx + bx; + v.y = t*my + by; } } // class LinearInterpolator private LinearMetrics() { } - public double distance(Point2D p1, Point2D p2) { + public double distance(Coordinate p1, Coordinate p2) { return p1.distance(p2); } - public Interpolator getInterpolator(Point2D p1, Point2D p2) { + public Interpolator getInterpolator(Coordinate p1, Coordinate p2) { return new LinearInterpolator(p1, p2); } } diff -r ee760729f6b8 -r aec85d00d82c gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearToMap.java --- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearToMap.java Tue Dec 15 15:58:58 2009 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearToMap.java Tue Dec 15 22:25:53 2009 +0000 @@ -1,6 +1,6 @@ package de.intevation.gnv.math; -import java.awt.geom.Point2D; +import com.vividsolutions.jts.geom.Coordinate; import java.util.List; import java.util.Iterator; @@ -11,20 +11,6 @@ */ public class LinearToMap { - public interface Interpolator { - - void interpolate(double t, Point2D v); - - } // interface Interpolator - - public interface Metrics { - - double distance(Point2D p1, Point2D p2); - - Interpolator getInterpolator(Point2D p1, Point2D p2); - - } // interface Metrics - public static final class Range { private Range next; @@ -32,8 +18,8 @@ private double to; private double b; - private Point2D p1; - private Point2D p2; + private Coordinate p1; + private Coordinate p2; private Interpolator interpolator; @@ -44,8 +30,8 @@ double from, double to, Interpolator interpolator, - Point2D p1, - Point2D p2 + Coordinate p1, + Coordinate p2 ) { this.from = from; this.to = to; @@ -58,7 +44,7 @@ : 1.0d/(to - from); } - public void eval(double x, Point2D v) { + public void eval(double x, Coordinate v) { interpolator.interpolate((x - from)*b, v); } @@ -66,11 +52,11 @@ return x >= from && x <= to; } - public Point2D startPoint() { + public Coordinate startPoint() { return p1; } - public Point2D endPoint() { + public Coordinate endPoint() { return p2; } } // class Range @@ -82,10 +68,10 @@ } public LinearToMap( - List path, - double from, - double to, - Metrics metrics + List path, + double from, + double to, + Metrics metrics ) { double diagramLength = Math.abs(to - from); @@ -96,8 +82,8 @@ Range last = null; for (int i = 1, N = path.size(); i < N; ++i) { - Point2D p1 = (Point2D)path.get(i-1); - Point2D p2 = (Point2D)path.get(i); + Coordinate p1 = path.get(i-1); + Coordinate p2 = path.get(i); double segmentLength = metrics.distance(p1, p2); double relativeLength = segmentLength / worldLength; @@ -139,7 +125,7 @@ return null; } - public boolean locate(double diagramX, Point2D v) { + public boolean locate(double diagramX, Coordinate v) { Range range = locateRange(diagramX); if (range == null) { return false; @@ -148,11 +134,14 @@ return true; } - public static double length(List path, Metrics metrics) { + public static double length( + List path, + Metrics metrics + ) { double sum = 0d; for (int i = path.size()-1; i >= 1; --i) { - Point2D p1 = (Point2D)path.get(i); - Point2D p2 = (Point2D)path.get(i-1); + Coordinate p1 = path.get(i); + Coordinate p2 = path.get(i-1); sum += metrics.distance(p1, p2); } return sum; diff -r ee760729f6b8 -r aec85d00d82c gnv-artifacts/src/main/java/de/intevation/gnv/math/Metrics.java --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Metrics.java Tue Dec 15 22:25:53 2009 +0000 @@ -0,0 +1,14 @@ +package de.intevation.gnv.math; + +import com.vividsolutions.jts.geom.Coordinate; + +/** + * @author Sascha L. Teichmann + */ +public interface Metrics +{ + double distance(Coordinate p1, Coordinate p2); + + Interpolator getInterpolator(Coordinate p1, Coordinate p2); +} +// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8: diff -r ee760729f6b8 -r aec85d00d82c gnv-artifacts/src/main/java/de/intevation/gnv/math/Point2d.java --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Point2d.java Tue Dec 15 22:25:53 2009 +0000 @@ -0,0 +1,95 @@ +package de.intevation.gnv.math; + +import java.util.Comparator; + +import com.vividsolutions.jts.geom.Envelope; +import com.vividsolutions.jts.geom.Coordinate; + +/** + * @author Sascha L. Teichmann + */ +public class Point2d +extends Coordinate +{ + public static final double EPSILON = 1e-5d; + + public static final Comparator X_COMPARATOR = new Comparator() { + public int compare(Object a, Object b) { + double xa = ((Coordinate)a).x; + double xb = ((Coordinate)b).x; + if (xa < xb) return -1; + if (xa > xb) return +1; + return 0; + } + }; + + public static final Comparator Y_COMPARATOR = new Comparator() { + public int compare(Object a, Object b) { + double ya = ((Coordinate)a).y; + double yb = ((Coordinate)b).y; + if (ya < yb) return -1; + if (ya > yb) return +1; + return 0; + } + }; + + public static class InverseL1Comparator + implements Comparator + { + private Point2d ref; + + public InverseL1Comparator(Point2d ref) { + this.ref = ref; + } + + public int compare(Object a, Object b) { + Point2d pa = (Point2d)a; + Point2d pb = (Point2d)b; + double da = ref.L1(pa); + double db = ref.L1(pb); + if (da < db) return -1; + if (da > db) return +1; + return 0; + } + } // class InverseL1Comparator + + public int i; + public int j; + + public Point2d() { + } + + public Point2d(double x, double y, double z, int i, int j) { + super(x, y, z); + this.i = i; + this.j = j; + } + + + public double L1(Point2d other) { + return L1(this, other); + } + + public static double L1(Coordinate a, Coordinate b) { + return Math.abs(a.x - b.x) + Math.abs(a.y - b.y); + } + + public Envelope envelope() { + return envelope(EPSILON); + } + + public Envelope envelope(double epsilon) { + return new Envelope( + x-epsilon, x+epsilon, + y-epsilon, y+epsilon); + } + + public boolean hasIGap(Point2d other) { + return Math.abs(i - other.i) > 1; + } + + public boolean hasJGap(Point2d other) { + return Math.abs(j - other.j) > 1; + } +} +// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8: diff -r ee760729f6b8 -r aec85d00d82c gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontal/HorizontalProfileMeshCrossOutputState.java --- a/gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontal/HorizontalProfileMeshCrossOutputState.java Tue Dec 15 15:58:58 2009 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontal/HorizontalProfileMeshCrossOutputState.java Tue Dec 15 22:25:53 2009 +0000 @@ -193,8 +193,10 @@ System.arraycopy(filterValues, 0, addedFilterValues, 0, filterValues.length); addedFilterValues[filterValues.length] = additionWhere; - result = queryExecutor.executeQuery(this.queryID, - addedFilterValues); + result = process( + queryExecutor.executeQuery( + this.queryID, + addedFilterValues)); } catch (ParseException e) { log.error(e,e); @@ -213,7 +215,12 @@ } return result; } - - + protected Collection process(Collection input) { + // TODO: split by additional parameters, interpolate the + // values, and create a new dummy result set. + + return input; + + } }