changeset 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 ee760729f6b8
children 1ab23cd66870
files gnv-artifacts/ChangeLog gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolator.java gnv-artifacts/src/main/java/de/intevation/gnv/math/L1Comparator.java gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearMetrics.java gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearToMap.java gnv-artifacts/src/main/java/de/intevation/gnv/math/Metrics.java gnv-artifacts/src/main/java/de/intevation/gnv/math/Point2d.java gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontal/HorizontalProfileMeshCrossOutputState.java
diffstat 9 files changed, 393 insertions(+), 51 deletions(-) [+]
line wrap: on
line diff
--- 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	<sascha.teichmann@intevation.de>
+
+	* 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	<sascha.teichmann@intevation.de>
 
 	* src/main/java/de/intevation/gnv/math/LinearToMap.java: Map linear
--- /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:
--- /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:
--- /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:
+
--- 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);
     }
 }
--- 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<? extends Coordinate> 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<? extends Coordinate> 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;
--- /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:
--- /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:
--- 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<Result> process(Collection<Result> input) {
+        // TODO: split by additional parameters, interpolate the
+        // values, and create a new dummy result set.
+
+        return input;
+
+    }
 }

http://dive4elements.wald.intevation.org