changeset 514:d9d933e06875

Fixed gnv/issue153 gnv-artifacts/trunk@608 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 22 Jan 2010 18:22:11 +0000
parents ca5048e4e515
children 234d9892e497
files gnv-artifacts/ChangeLog gnv-artifacts/src/main/java/de/intevation/gnv/chart/TimeSeriesChart.java gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java gnv-artifacts/src/main/java/de/intevation/gnv/math/GridCell.java gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearFunction.java gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontal/HorizontalProfileOutputState.java
diffstat 7 files changed, 262 insertions(+), 124 deletions(-) [+]
line wrap: on
line diff
--- a/gnv-artifacts/ChangeLog	Fri Jan 22 15:45:59 2010 +0000
+++ b/gnv-artifacts/ChangeLog	Fri Jan 22 18:22:11 2010 +0000
@@ -1,3 +1,33 @@
+2010-01-22	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
+
+	* src/main/java/de/intevation/gnv/math/GridCell.java: New.
+	  A 4-tupel of neighbored points in the mesh. It is valid
+	  to interpolate in this area.
+
+	* src/main/java/de/intevation/gnv/math/AreaInterpolation.java:
+	  The algorithm how neighbored points in the mesh are determined
+	  has changed. Now all incoming points are tiled into GridCells.
+	  If there are gaps in i,j the corresponding tile is omited.
+	  These tiles are stored in an R tree. To lookup a point in
+	  world coordinates the spatial index is queried. If no result
+	  is found the point is ignore as a gap. If a fitting grid cell
+	  is found the interpolation in done between the four points
+	  of that cell is performed. Special gap checking is not needed any 
+	  longer. This fixes gnv/issue153 because there are no assumptions
+	  about axis aligned points any more.
+
+	* src/main/java/de/intevation/gnv/math/Interpolation2D.java: Used
+	  euclid distance to estimate spatial buffer size. TODO: Remove
+	  this code when adjusting the "Profilschnitte" to the same logic
+	  as used in "Horzontalschnitte" now.
+
+	* src/main/java/de/intevation/gnv/math/LinearFunction.java: Added
+	  author's email.
+
+	* src/main/java/de/intevation/gnv/chart/TimeSeriesChart.java,
+	  src/main/java/de/intevation/gnv/state/profile/horizontal/HorizontalProfileOutputState.java: 
+	  Cleanup imports.
+
 2010-01-22  Ingo Weinzierl <ingo.weinzierl@intevation.de>
 
 	* src/main/java/de/intevation/gnv/artifacts/GNVArtifactBase.java: Removed
@@ -27,6 +57,7 @@
 	  Added the possibility to configure the different width of the Meshes.
 
 2010-01-22  Tim Englich  <tim.englich@intevation.de>
+
 	* src/test/java/de/intevation/gnv/artifacts/TestArtifactDatabase.java (serviceNamesAndDescriptions),(process):
 	  Fixed Compiler-Error after adding further Methods to the Interface ArtifactDatabase.
 
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/chart/TimeSeriesChart.java	Fri Jan 22 15:45:59 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/chart/TimeSeriesChart.java	Fri Jan 22 18:22:11 2010 +0000
@@ -1,7 +1,13 @@
 package de.intevation.gnv.chart;
 
-import java.text.DateFormat;
-import java.text.SimpleDateFormat;
+import de.intevation.gnv.artifacts.ressource.RessourceFactory;
+
+import de.intevation.gnv.geobackend.base.Result;
+
+import de.intevation.gnv.state.describedata.KeyValueDescibeData;
+
+import de.intevation.gnv.timeseries.gap.TimeGap;
+
 import java.util.Collection;
 import java.util.Date;
 import java.util.HashMap;
@@ -11,25 +17,23 @@
 
 import org.apache.log4j.Logger;
 
+import org.jfree.chart.ChartFactory;
 import org.jfree.chart.ChartTheme;
-import org.jfree.chart.ChartFactory;
+
 import org.jfree.chart.axis.Axis;
 import org.jfree.chart.axis.DateAxis;
-import org.jfree.chart.plot.XYPlot;
+
 import org.jfree.chart.plot.PlotOrientation;
+import org.jfree.chart.plot.XYPlot;
+
 import org.jfree.data.general.Series;
+
+import org.jfree.data.time.Minute;
 import org.jfree.data.time.TimeSeries;
-import org.jfree.data.time.Minute;
 import org.jfree.data.time.TimeSeriesCollection;
 
-import de.intevation.gnv.artifacts.ressource.RessourceFactory;
-import de.intevation.gnv.geobackend.base.Result;
-import de.intevation.gnv.state.describedata.KeyValueDescibeData;
-import de.intevation.gnv.timeseries.gap.TimeGap;
-
-
 /**
- * @author Ingo Weinzierl <ingo.weinzierl@intevation.de>
+ * @author Ingo Weinzierl (ingo.weinzierl@intevation.de)
  */
 public class TimeSeriesChart
 extends      AbstractXYLineChart
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java	Fri Jan 22 15:45:59 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java	Fri Jan 22 18:22:11 2010 +0000
@@ -3,15 +3,13 @@
 import com.vividsolutions.jts.geom.Coordinate;
 import com.vividsolutions.jts.geom.Envelope;
 
-import com.vividsolutions.jts.index.quadtree.Quadtree;
+import com.vividsolutions.jts.index.strtree.STRtree;
 
 import java.awt.Dimension;
 
 import java.io.Serializable;
 
-import java.util.ArrayList;
 import java.util.Arrays;
-import java.util.Collections;
 import java.util.List;
 
 import org.apache.log4j.Logger;
@@ -44,6 +42,9 @@
         return raster;
     }
 
+    public static final double EPS = 1e-6d;
+
+
     public boolean interpolate(
         List<? extends Point2d> points,
         Envelope                boundingBox,
@@ -53,9 +54,14 @@
         boolean debug = log.isDebugEnabled();
 
         if (points == null || points.isEmpty()) {
-            if (debug) {
-                log.debug("no points to interpolate");
-            }
+            log.warn("no points to interpolate");
+            return false;
+        }
+
+        List<GridCell> cells = GridCell.pointsToGridCells(points);
+
+        if (cells.isEmpty()) {
+            log.warn("no cells to interpolate");
             return false;
         }
 
@@ -64,47 +70,25 @@
 
         double cellWidth  = boundingBox.getWidth()  / W;
         double cellHeight = boundingBox.getHeight() / H;
-        Envelope gridEnvelope = new Envelope(boundingBox);
-
-        gridEnvelope.expandBy(2d*cellWidth, 2d*cellWidth);
-
-        ArrayList<Point2d>relevantPoints =
-            new ArrayList<Point2d>();
 
-        for (int i = points.size()-1; i >= 0; --i) {
-            if (gridEnvelope.contains(points.get(i))) {
-                relevantPoints.add(points.get(i));
-            }
-        }
+        STRtree spatialIndex = new STRtree();
 
-        double [] buffer = Interpolation2D
-            .calculateBufferRemoveOutlier(relevantPoints);
-        double dxMax = buffer[0];
-        double dyMax = buffer[1];
-
-        Quadtree spatialIndex = new Quadtree();
-
-        for (int i = relevantPoints.size()-1; i >= 0; --i) {
-            Point2d p = relevantPoints.get(i);
-            spatialIndex.insert(p.envelope(), p);
+        for (GridCell cell: cells) {
+            spatialIndex.insert(cell.getEnvelope(), cell);
         }
 
         if (debug) {
-            log.debug("points: "        + relevantPoints.size());
             log.debug("width:  "        + boundingBox.getWidth());
             log.debug("height:  "       + boundingBox.getHeight());
             log.debug("sample width:  " + W);
             log.debug("sample height: " + H);
             log.debug("cell width:  "   + cellWidth);
             log.debug("cell height: "   + cellHeight);
-            log.debug("buffer x: "      + dxMax);
-            log.debug("buffer y: "      + dyMax);
         }
 
-        Envelope     queryBuffer = new Envelope();
-        Point2d []   neighbors   = new Point2d[4];
-        Coordinate   center      = new Coordinate();
-        L1Comparator invL1       = new L1Comparator(center);
+        Envelope   queryBuffer     = new Envelope();
+        Coordinate center          = new Coordinate();
+        GridCell.CellFinder finder = new GridCell.CellFinder();
 
         double [] raster = new double[W*H];
         Arrays.fill(raster, Double.NaN);
@@ -119,64 +103,38 @@
             double y = j*cellHeight + 0.5d*cellHeight + minY;
             double x = 0.5d*cellWidth + minX;
             for (int i = 0; i < W; ++i, x += cellWidth) {
-                queryBuffer.init(
-                    x - dxMax, x + dxMax, 
-                    y - dyMax, y + dyMax);
 
-                List potential = spatialIndex.query(queryBuffer);
+                queryBuffer.init(x - EPS, x + EPS, y - EPS, y + EPS);
+                center.x = x; center.y = y;
+                finder.prepare(center);
+                spatialIndex.query(queryBuffer, finder);
 
-                if (potential.size() < 4) {
+                GridCell found = finder.found;
+
+                if (found == null) {
                     continue;
                 }
 
-                center.x = x; center.y = y;
-                Collections.sort(potential, invL1);
-
-                neighbors[0] = neighbors[1] = 
-                neighbors[2] = neighbors[3] = null;
-
-                int mask = 0;
-                // reversed order is essential here!
-                for (int k = potential.size()-1; k >= 0; --k) {
-                    Point2d n = (Point2d)potential.get(k);
-                    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
-                // we do not have any gaps and we are below surface.
-                if (numNeighbors == 4
-                && !neighbors[0].hasIGap(neighbors[1])
-                && !neighbors[1].hasJGap(neighbors[3])
-                && !neighbors[3].hasIGap(neighbors[2])
-                && !neighbors[2].hasJGap(neighbors[0])
-                && depth.depth(center) <= 0d
-                ) {
-                    double z1 = Interpolation2D.interpolate(
-                        neighbors[0].x, neighbors[0].z,
-                        neighbors[1].x, neighbors[1].z,
-                        center.x);
-                    double z2 = Interpolation2D.interpolate(
-                        neighbors[2].x, neighbors[2].z,
-                        neighbors[3].x, neighbors[3].z,
-                        center.x);
-                    double y1 = Interpolation2D.interpolate(
-                        neighbors[0].x, neighbors[0].y,
-                        neighbors[1].x, neighbors[1].y,
-                        center.x);
-                    double y2 = Interpolation2D.interpolate(
-                        neighbors[2].x, neighbors[2].y,
-                        neighbors[3].x, neighbors[3].y,
-                        center.x);
-                    raster[row+i] = Interpolation2D.interpolate(
-                        y1, z1,
-                        y2, z2,
-                        center.y);
-                }
+                double z1 = Interpolation2D.interpolate(
+                    found.p1.x, found.p1.z,
+                    found.p2.x, found.p2.z,
+                    center.x);
+                double z2 = Interpolation2D.interpolate(
+                    found.p3.x, found.p3.z,
+                    found.p4.x, found.p4.z,
+                    center.x);
+                double y1 = Interpolation2D.interpolate(
+                    found.p1.x, found.p1.y,
+                    found.p2.x, found.p2.y,
+                    center.x);
+                double y2 = Interpolation2D.interpolate(
+                    found.p3.x, found.p3.y,
+                    found.p4.x, found.p4.y,
+                    center.x);
+                raster[row+i] = Interpolation2D.interpolate(
+                    y1, z1,
+                    y2, z2,
+                    center.y);
             }
         }
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/GridCell.java	Fri Jan 22 18:22:11 2010 +0000
@@ -0,0 +1,134 @@
+package de.intevation.gnv.math;
+
+import com.vividsolutions.jts.geom.Coordinate;
+import com.vividsolutions.jts.geom.Envelope;
+import com.vividsolutions.jts.geom.Geometry;
+import com.vividsolutions.jts.geom.GeometryFactory;
+import com.vividsolutions.jts.geom.LinearRing;
+import com.vividsolutions.jts.geom.Point;
+import com.vividsolutions.jts.geom.Polygon;
+
+import com.vividsolutions.jts.index.ItemVisitor;
+
+import java.io.Serializable;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+
+/**
+ *  @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
+ */
+public class GridCell
+implements   Serializable
+{
+    public static final class CellFinder
+    implements                ItemVisitor
+    {
+        public    GridCell found;
+
+        protected Point    point;
+
+        public CellFinder() {
+        }
+
+        public void prepare(Coordinate center) {
+            found = null;
+            point = GEOMETRY_FACTORY.createPoint(center);
+        }
+
+        public void visitItem(Object item) {
+            if (found == null) {
+                GridCell cell = (GridCell)item;
+                if (cell.contains(point)) {
+                    found = cell;
+                }
+            }
+        }
+    } // class CellFinder
+
+    public Point2d p1;
+    public Point2d p2;
+    public Point2d p3;
+    public Point2d p4;
+
+    protected Polygon polygon;
+
+    public static final GeometryFactory GEOMETRY_FACTORY
+        = new GeometryFactory();
+
+    public GridCell() {
+    }
+
+    public GridCell(Point2d p1, Point2d p2, Point2d p3, Point2d p4) {
+        this.p1 = p1;
+        this.p2 = p2;
+        this.p3 = p3;
+        this.p4 = p4;
+        createPolygon();
+    }
+
+    protected void createPolygon() {
+        LinearRing shell = GEOMETRY_FACTORY
+            .createLinearRing(new Coordinate [] { p1, p2, p3, p4, p1 });
+        polygon = GEOMETRY_FACTORY.createPolygon(shell, null);
+    }
+
+    public Envelope getEnvelope() {
+        return polygon.getEnvelopeInternal();
+    }
+
+    public boolean contains(Geometry coord) {
+        return polygon.contains(coord);
+    }
+
+    public static List<GridCell> pointsToGridCells(
+        List<? extends Point2d> points
+    ) {
+        int minI = Integer.MAX_VALUE;
+        int maxI = Integer.MIN_VALUE;
+        int minJ = Integer.MAX_VALUE;
+        int maxJ = Integer.MIN_VALUE;
+
+        HashMap<Integer, HashMap<Integer, Point2d>> rows =
+            new HashMap<Integer, HashMap<Integer, Point2d>>();
+
+        for (Point2d p: points) {
+
+            if (p.i < minI) minI = p.i;
+            if (p.i > maxI) maxI = p.i;
+            if (p.j < minJ) minJ = p.j;
+            if (p.j > maxJ) maxJ = p.j;
+
+            HashMap<Integer, Point2d> row = rows.get(p.i);
+
+            if (row == null) {
+                rows.put(p.i, row = new HashMap<Integer, Point2d>());
+            }
+
+            row.put(p.j, p);
+        }
+
+        ArrayList<GridCell> cells = new ArrayList<GridCell>(points.size());
+
+        for (int i = minI + 1; i <= maxI; ++i) {
+            HashMap<Integer, Point2d> row1 = rows.get(i-1);
+            HashMap<Integer, Point2d> row2 = rows.get(i);
+            if (row1 != null && row2 != null) {
+                for (int j = minJ + 1; j < maxJ; ++j) {
+                    Point2d p1 = row1.get(j-1);
+                    Point2d p2 = row1.get(j);
+                    Point2d p3 = row2.get(j);
+                    Point2d p4 = row2.get(j-1);
+
+                    if (p1 != null && p2 != null &&  p3 != null && p4 != null) {
+                        cells.add(new GridCell(p1, p2, p3, p4));
+                    }
+                } // for all columns in row
+            }
+        } // for all rows
+
+        return cells;
+    }
+}
+// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Fri Jan 22 15:45:59 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Fri Jan 22 18:22:11 2010 +0000
@@ -72,7 +72,8 @@
         for (ArrayList<Point2d> v: jMap.values()) {
             Collections.sort(v, Point2d.Y_COMPARATOR);
             for (int i = 1, L = v.size(); i < L; ++i) {
-                double dy = Math.abs(v.get(i).y - v.get(i-1).y);
+                // double dy = Math.abs(v.get(i).y - v.get(i-1).y);
+                double dy = v.get(i).distance(v.get(i-1));
                 deltas.add(dy);
                 mean.increment(dy);
                 sd  .increment(dy);
@@ -104,7 +105,8 @@
         for (ArrayList<Point2d> v: iMap.values()) {
             Collections.sort(v, Point2d.X_COMPARATOR);
             for (int i = 1, L = v.size(); i < L; ++i) {
-                double dx = Math.abs(v.get(i).x - v.get(i-1).x);
+                //double dx = Math.abs(v.get(i).x - v.get(i-1).x);
+                double dx = v.get(i).distance(v.get(i-1));
                 deltas.add(dx);
                 mean.increment(dx);
                 sd  .increment(dx);
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearFunction.java	Fri Jan 22 15:45:59 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/LinearFunction.java	Fri Jan 22 18:22:11 2010 +0000
@@ -7,7 +7,7 @@
 import org.apache.commons.math.analysis.UnivariateRealFunction;
 
 /**
- *  @author Sascha L. Teichmann
+ *  @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
  */
 public class LinearFunction
 implements   ParametricRealFunction
@@ -54,4 +54,4 @@
         return new double [] { x, 1f };
     }
 }
-// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8:
+// 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/HorizontalProfileOutputState.java	Fri Jan 22 15:45:59 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontal/HorizontalProfileOutputState.java	Fri Jan 22 18:22:11 2010 +0000
@@ -3,14 +3,41 @@
  */
 package de.intevation.gnv.state.profile.horizontal;
 
-import com.vividsolutions.jts.geom.Point;
+import com.vividsolutions.jts.io.ParseException;
 import com.vividsolutions.jts.io.WKTReader;
-import com.vividsolutions.jts.io.ParseException;
+
+import de.intevation.artifacts.CallContext;
+
+import de.intevation.gnv.chart.Chart;
+import de.intevation.gnv.chart.ChartLabels;
+import de.intevation.gnv.chart.HorizontalProfileChart;
+
+import de.intevation.gnv.exports.DefaultExport;
+import de.intevation.gnv.exports.DefaultProfile;
+
+import de.intevation.gnv.exports.Export.Profile;
+
+import de.intevation.gnv.exports.ShapeDataCollector;
+
+import de.intevation.gnv.geobackend.base.Result;
+
+import de.intevation.gnv.state.describedata.KeyValueDescibeData;
+
+import de.intevation.gnv.state.exception.StateException;
+
+import de.intevation.gnv.state.timeseries.TimeSeriesOutputState;
+
+import de.intevation.gnv.statistics.HorizontalProfileStatistics;
+import de.intevation.gnv.statistics.Statistics;
+
+import de.intevation.gnv.utils.WKTUtils;
 
 import java.io.IOException;
 import java.io.OutputStream;
 import java.io.UnsupportedEncodingException;
+
 import java.text.SimpleDateFormat;
+
 import java.util.Collection;
 import java.util.Date;
 import java.util.Iterator;
@@ -20,26 +47,8 @@
 
 import org.jfree.chart.ChartTheme;
 
-import de.intevation.artifacts.CallContext;
-
-import de.intevation.gnv.chart.Chart;
-import de.intevation.gnv.chart.ChartLabels;
-import de.intevation.gnv.chart.HorizontalProfileChart;
-import de.intevation.gnv.exports.DefaultExport;
-import de.intevation.gnv.exports.ShapeDataCollector;
-import de.intevation.gnv.exports.DefaultProfile;
-import de.intevation.gnv.exports.Export.Profile;
-import de.intevation.gnv.geobackend.base.Result;
-import de.intevation.gnv.state.describedata.KeyValueDescibeData;
-import de.intevation.gnv.state.exception.StateException;
-import de.intevation.gnv.state.timeseries.TimeSeriesOutputState;
-import de.intevation.gnv.statistics.HorizontalProfileStatistics;
-import de.intevation.gnv.statistics.Statistics;
-import de.intevation.gnv.utils.WKTUtils;
-
 /**
- * @author Tim Englich <tim.englich@intevation.de>
- * 
+ * @author Tim Englich (tim.englich@intevation.de)
  */
 public class HorizontalProfileOutputState
 extends      TimeSeriesOutputState

http://dive4elements.wald.intevation.org