# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1264281405 0
# Node ID 4e347624ee7cd276780a8f9b525bdabdfaf10477
# Parent  464e03bf786baa2b2abf47e476b1351bd9bc4647
Last part to fix gnv/issue153. Now 'Profilschnitte', 'Horizontalschnitte' and 'horizontale Schnittprofile'
all use the same x/y interpolation code.



gnv-artifacts/trunk@613 c6561f87-3c4e-4783-a992-168aeb5c3f6f

diff -r 464e03bf786b -r 4e347624ee7c gnv-artifacts/ChangeLog
--- a/gnv-artifacts/ChangeLog	Sat Jan 23 18:10:34 2010 +0000
+++ b/gnv-artifacts/ChangeLog	Sat Jan 23 21:16:45 2010 +0000
@@ -1,3 +1,27 @@
+2010-01-23	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
+
+	* src/main/java/de/intevation/gnv/state/profile/horizontal/HorizontalProfileMeshCrossOutputState.java:
+	  Fixed bug when accessing i and j columns of SQL dataset. This
+	  prevented gap detection in "horizontale Schnittprofile" from working.
+
+	* src/main/java/de/intevation/gnv/math/Interpolation2D.java:
+	  "horizontale Schnittprofile" are now using the grid cell mechanism
+	  too. This should fix all remaining problems to solve gnv/issue153.
+	  The culling of too much points is controlled with the system property
+	  "gnv.interpolation2d.cull.point.threshold" with the same semantics 
+	  as in 'Profilschnitt' and 'Horizontalschnitt'.
+	  The spatial buffer size estimation code is removed because it is 
+	  not needed any longer.
+
+	* src/main/java/de/intevation/gnv/math/Interpolation3D.java: Moved some
+	  code to Interpolation2D.
+
+	* src/main/java/de/intevation/gnv/math/GridCell.java: Added some
+	  debug information about the number of found cells.
+
+	* src/main/java/de/intevation/gnv/utils/WKTUtils.java:
+	  Cleanup imports.
+
 2010-01-23	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
 
 	* src/main/java/de/intevation/gnv/math/Interpolation3D.java:
@@ -27,7 +51,7 @@
 
 2010-01-23	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
 
-	* contrib/palette2qgis.xsl: New. XST tranformation to turn a
+	* contrib/palette2qgis.xsl: New. XST transformation to turn a
 	  palette XML file into a style definition suitable to be used
 	  in QGIS. Tested with QGIS 1.4.0-Enceladus. Usage:
 
diff -r 464e03bf786b -r 4e347624ee7c gnv-artifacts/src/main/java/de/intevation/gnv/math/GridCell.java
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/GridCell.java	Sat Jan 23 18:10:34 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/GridCell.java	Sat Jan 23 21:16:45 2010 +0000
@@ -127,10 +127,6 @@
             row.put(p.j, p);
         }
 
-        if (log.isDebugEnabled()) {
-            log.debug("culled points: " + culled);
-        }
-
         ArrayList<GridCell> cells = new ArrayList<GridCell>(points.size());
 
         for (int i = minI + 1; i <= maxI; ++i) {
@@ -150,6 +146,13 @@
             }
         } // for all rows
 
+        if (log.isDebugEnabled()) {
+            log.debug("culled points: " + culled);
+            log.debug("min/max i: " + minI + " / " + maxI);
+            log.debug("min/max j: " + minJ + " / " + maxJ);
+            log.debug("cells found: " + cells.size());
+        }
+
         return cells;
     }
 }
diff -r 464e03bf786b -r 4e347624ee7c gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Sat Jan 23 18:10:34 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Sat Jan 23 21:16:45 2010 +0000
@@ -3,18 +3,10 @@
 import com.vividsolutions.jts.geom.Coordinate;
 import com.vividsolutions.jts.geom.Envelope;
 
-import com.vividsolutions.jts.index.quadtree.Quadtree;
-
-import gnu.trove.TDoubleArrayList;
+import com.vividsolutions.jts.index.strtree.STRtree;
 
-import java.util.ArrayList;
-import java.util.Collections;
-import java.util.HashMap;
 import java.util.List;
 
-import org.apache.commons.math.stat.descriptive.moment.Mean;
-import org.apache.commons.math.stat.descriptive.moment.StandardDeviation;
-
 import org.apache.log4j.Logger;
 
 /**
@@ -24,6 +16,11 @@
 {
     private static Logger log = Logger.getLogger(Interpolation2D.class);
 
+    public static final int CULL_POINT_THRESHOLD = Integer.getInteger(
+        "gnv.interpolation2d.cull.point.threshold", 1000);
+
+    public static final double EPS = 1e-6d;
+
     public interface Consumer {
         void interpolated(Coordinate point, boolean success);
     } // interface Consumer
@@ -31,164 +28,22 @@
     private Interpolation2D() {
     }
 
-    private static final double dist(double a, double b) {
-        a -= b;
-        return Math.sqrt(a*a);
-    }
-
-    private static final double OUTLIER_EPS = 0.4d;
-
-    public static final double [] calculateBufferRemoveOutlier(
-        List <? extends Point2d> points
+    public static Envelope pathBoundingBox(
+        List<? extends Coordinate> path,
+        double extra
     ) {
-        HashMap<Integer, ArrayList<Point2d>> iMap = 
-            new HashMap<Integer, ArrayList<Point2d>>();
-
-        HashMap<Integer, ArrayList<Point2d>> jMap = 
-            new HashMap<Integer, ArrayList<Point2d>>();
-
-        for (int k = points.size()-1; k >= 0; --k) {
-            Point2d p = points.get(k);
-
-            ArrayList<Point2d> jList = jMap.get(p.j);
-            ArrayList<Point2d> iList = iMap.get(p.i);
-
-            if (jList == null) {
-                jMap.put(p.j, jList = new ArrayList<Point2d>());
-            }
-            jList.add(p);
+        int N = path.size();
+        Envelope area = new Envelope(path.get(N-1));
 
-            if (iList == null) {
-                iMap.put(p.i, iList = new ArrayList<Point2d>());
-            }
-            iList.add(p);
-        }
-
-        TDoubleArrayList deltas = new TDoubleArrayList();
-
-        Mean mean            = new Mean();
-        StandardDeviation sd = new StandardDeviation();
-
-        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 = v.get(i).distance(v.get(i-1));
-                deltas.add(dy);
-                mean.increment(dy);
-                sd  .increment(dy);
-            }
-        }
-
-        deltas.sort();
-
-        double m   = mean.getResult();
-        double s   = sd  .getResult();
-        double eps = OUTLIER_EPS*s;
-
-        int yi = deltas.size() - 1;
-        while (yi > 0 && dist(deltas.get(yi), m) > eps) {
-            --yi;
-        }
-
-        double dyMax = deltas.get(yi) + 1e-5d;
-
-        if (log.isDebugEnabled()) {
-            log.debug("mean  y: " + m);
-            log.debug("sigma y: " + s);
+        for (int i = N-2; i >= 0; --i) {
+            area.expandToInclude(path.get(i));
         }
 
-        deltas.reset();
-        mean.clear();
-        sd  .clear();
-
-        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 = v.get(i).distance(v.get(i-1));
-                deltas.add(dx);
-                mean.increment(dx);
-                sd  .increment(dx);
-            }
-        }
-
-        deltas.sort();
-
-        m   = mean.getResult();
-        s   = sd  .getResult();
-        eps = OUTLIER_EPS*s;
-
-        int xi = deltas.size() - 1;
-
-        while (xi > 0 && dist(deltas.get(xi), m) > eps) {
-            --xi;
-        }
-
-        double dxMax = deltas.get(xi) + 1e-5d;
-
-        if (log.isDebugEnabled()) {
-            log.debug("mean  y: " + m);
-            log.debug("sigma y: " + s);
-        }
-
-        return new double [] { dxMax, dyMax };
-    }
-
-    public static final double [] calculateBuffer(
-        List <? extends Point2d> points
-    ) {
-        HashMap<Integer, ArrayList<Point2d>> iMap = 
-            new HashMap<Integer, ArrayList<Point2d>>();
+        area.expandBy(
+            extra*area.getWidth(),
+            extra*area.getHeight());
 
-        HashMap<Integer, ArrayList<Point2d>> jMap = 
-            new HashMap<Integer, ArrayList<Point2d>>();
-
-        for (int k = points.size()-1; k >= 0; --k) {
-            Point2d p = points.get(k);
-
-            ArrayList<Point2d> jList = jMap.get(p.j);
-            ArrayList<Point2d> iList = iMap.get(p.i);
-
-            if (jList == null) {
-                jMap.put(p.j, jList = new ArrayList<Point2d>());
-            }
-            jList.add(p);
-
-            if (iList == null) {
-                iMap.put(p.i, iList = new ArrayList<Point2d>());
-            }
-            iList.add(p);
-        }
-
-        double dxMax = -Double.MAX_VALUE;
-        double dyMax = -Double.MAX_VALUE;
-
-        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);
-                if (dy > dyMax) {
-                    dyMax = dy;
-                }
-            }
-        }
-
-        dyMax += 1e-5d;
-
-        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);
-                if (dx > dxMax) {
-                    dxMax = dx;
-                }
-            }
-        }
-
-        dxMax += 1e-5d;
-
-        return new double [] { dxMax, dyMax };
+        return area;
     }
 
     public static void interpolate(
@@ -214,20 +69,23 @@
             return;
         }
 
-        double [] buffer = calculateBuffer(points);
-        double dxMax = buffer[0];
-        double dyMax = buffer[1];
+        Envelope relevantArea = M > CULL_POINT_THRESHOLD
+			? pathBoundingBox(path, 0.05d)
+			: null;
+        
+        List<GridCell> cells = GridCell.pointsToGridCells(
+            points, relevantArea);
 
-        if (debug) {
-            log.debug("buffer size: " + dxMax + " / " + dyMax);
+        if (cells.isEmpty()) {
+            log.warn("no cells found");
+            return;
         }
 
         // put into spatial index to speed up finding neighbors.
-        Quadtree spatialIndex = new Quadtree();
+        STRtree spatialIndex = new STRtree();
 
-        for (int i = 0; i < M; ++i) {
-            Point2d p = points.get(i);
-            spatialIndex.insert(p.envelope(), p);
+        for (GridCell cell: cells) {
+            spatialIndex.insert(cell.getEnvelope(), cell);
         }
 
         LinearToMap linearToMap = new LinearToMap(
@@ -235,85 +93,59 @@
 
         double dP = (to - from)/steps;
 
-        Coordinate center = new Coordinate();
-
-        Envelope queryBuffer = new Envelope();
-
-        Point2d [] neighbors = new Point2d[4];
+        Coordinate          center      = new Coordinate();
+        Envelope            queryBuffer = new Envelope();
+        GridCell.CellFinder finder      = new GridCell.CellFinder();
 
         int missedInterpolations = 0;
-        int interpolations = 0;
-
-        L1Comparator invL1 = new L1Comparator(center);
+        int interpolations       = 0;
 
         for (double p = from; p <= to; 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);
-
-            Collections.sort(potential, invL1);
-
-            neighbors[0] = neighbors[1] = 
-            neighbors[2] = neighbors[3] = null;
+                center.x - EPS, center.x + EPS, 
+                center.y - EPS, center.y + EPS);
 
-            /* bit code of neighbors
-               0---1
-               | x |
-               2---3
-            */
+            finder.prepare(center);
+            spatialIndex.query(queryBuffer, finder);
 
-            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;
+            GridCell found = finder.found;
+
+            if (found == null) {
+                consumer.interpolated(center, false);
+                ++missedInterpolations;
+                continue;
             }
 
-            int numNeighbors = Integer.bitCount(mask);
+            Point2d n0 = found.p1;
+            Point2d n1 = found.p2;
+            Point2d n2 = found.p3;
+            Point2d n3 = found.p4;
 
-            // only interpolate if we have all four neighbors
-            // and we do not have any gaps.
-            if (numNeighbors == 4
-            && !neighbors[0].hasIGap(neighbors[1])
-            && !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, true);
-                ++interpolations;
-            }
-            else {
-                consumer.interpolated(center, false);
-                ++missedInterpolations;
-            }
+            double z1 = interpolate(
+                n0.x, n0.z,
+                n1.x, n1.z,
+                center.x);
+            double z2 = interpolate(
+                n2.x, n2.z,
+                n3.x, n3.z,
+                center.x);
+            double y1 = interpolate(
+                n0.x, n0.y,
+                n1.x, n1.y,
+                center.x);
+            double y2 = interpolate(
+                n2.x, n2.y,
+                n3.x, n3.y,
+                center.x);
+            center.z = interpolate(
+                y1, z1,
+                y2, z2,
+                center.y);
+            consumer.interpolated(center, true);
+            ++interpolations;
         }
 
         if (debug) {
diff -r 464e03bf786b -r 4e347624ee7c gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java	Sat Jan 23 18:10:34 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java	Sat Jan 23 21:16:45 2010 +0000
@@ -98,21 +98,10 @@
             return false;
         }
 
-        Envelope relevantArea = null;
+        Envelope relevantArea = M > CULL_POINT_THRESHOLD
+			? Interpolation2D.pathBoundingBox(path, 0.05d)
+			: null;
         
-        if (M > CULL_POINT_THRESHOLD) {
-
-			relevantArea = new Envelope(path.get(N-1));
-
-			for (int i = N-2; i >= 0; --i) {
-				relevantArea.expandToInclude(path.get(i));
-			}
-
-            relevantArea.expandBy(
-                0.05d*relevantArea.getWidth(),
-                0.05d*relevantArea.getHeight());
-        }
-
         List<GridCell> cells = GridCell.pointsToGridCells(
             points, relevantArea);
 
diff -r 464e03bf786b -r 4e347624ee7c 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	Sat Jan 23 18:10:34 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontal/HorizontalProfileMeshCrossOutputState.java	Sat Jan 23 21:16:45 2010 +0000
@@ -3,46 +3,51 @@
  */
 package de.intevation.gnv.state.profile.horizontal;
 
-import java.util.Arrays;
-import java.util.ArrayList;
-import java.util.Collection;
-import java.util.List;
-import java.util.Locale;
-
-import org.apache.log4j.Logger;
-import org.w3c.dom.Node;
+import com.vividsolutions.jts.geom.Coordinate;
 
 import de.intevation.artifactdatabase.Config;
 
+import de.intevation.artifacts.CallContext;
+
 import de.intevation.gnv.artifacts.cache.CacheFactory;
 
+import de.intevation.gnv.artifacts.context.GNVArtifactContext;
+
+import de.intevation.gnv.chart.Chart;
+import de.intevation.gnv.chart.ChartLabels;
+import de.intevation.gnv.chart.HorizontalCrossProfileChart;
+
 import de.intevation.gnv.geobackend.base.DefaultResult;
 import de.intevation.gnv.geobackend.base.DefaultResultDescriptor;
 import de.intevation.gnv.geobackend.base.Result;
 import de.intevation.gnv.geobackend.base.ResultDescriptor;
+
 import de.intevation.gnv.geobackend.base.query.QueryExecutor;
 import de.intevation.gnv.geobackend.base.query.QueryExecutorFactory;
+
 import de.intevation.gnv.geobackend.base.query.exception.QueryException;
+
 import de.intevation.gnv.math.Interpolation2D;
 import de.intevation.gnv.math.LinearMetrics;
 import de.intevation.gnv.math.Point2d;
-import de.intevation.gnv.chart.Chart;
-import de.intevation.gnv.chart.ChartLabels;
-import de.intevation.gnv.chart.HorizontalCrossProfileChart;
 
 import de.intevation.gnv.state.InputData;
 
 import de.intevation.gnv.utils.DistanceCalculator;
+import de.intevation.gnv.utils.StringUtils;
 import de.intevation.gnv.utils.WKTUtils;
-import de.intevation.gnv.utils.StringUtils;
 
-import de.intevation.gnv.artifacts.context.GNVArtifactContext;
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.Collection;
+import java.util.List;
+import java.util.Locale;
 
-import de.intevation.artifacts.CallContext;
+import org.apache.log4j.Logger;
 
 import org.jfree.chart.ChartTheme;
 
-import com.vividsolutions.jts.geom.Coordinate;
+import org.w3c.dom.Node;
 
 /**
  * @author Tim Englich         (tim.englich@intevation.de)
@@ -360,8 +365,8 @@
             Coordinate coordinate =
                 WKTUtils.toCoordinate(result.getString("SHAPE"));
             double value = result.getDouble("YORDINATE");
-            int iPos     = result.getInteger("MEDIAN.MESHPOINT.JPOSITION");
-            int jPos     = result.getInteger("MEDIAN.MESHPOINT.JPOSITION");
+            int iPos     = result.getInteger("IPOSITION");
+            int jPos     = result.getInteger("JPOSITION");
             Point2d p = new Point2d(
                 coordinate.x,
                 coordinate.y,
diff -r 464e03bf786b -r 4e347624ee7c gnv-artifacts/src/main/java/de/intevation/gnv/utils/WKTUtils.java
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/utils/WKTUtils.java	Sat Jan 23 18:10:34 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/utils/WKTUtils.java	Sat Jan 23 21:16:45 2010 +0000
@@ -1,32 +1,39 @@
 package de.intevation.gnv.utils;
 
 import com.vividsolutions.jts.geom.Coordinate;
+import com.vividsolutions.jts.geom.LineString;
 import com.vividsolutions.jts.geom.Point;
 import com.vividsolutions.jts.geom.Polygon;
-import com.vividsolutions.jts.geom.LineString;
+
 import com.vividsolutions.jts.io.ParseException;
-
 import com.vividsolutions.jts.io.WKTReader;
 
+import de.intevation.gnv.artifacts.ressource.RessourceFactory;
+
 import de.intevation.gnv.geobackend.base.Result;
+
 import de.intevation.gnv.geobackend.base.query.QueryExecutor;
 import de.intevation.gnv.geobackend.base.query.QueryExecutorFactory;
+
 import de.intevation.gnv.geobackend.base.query.exception.QueryException;
-import de.intevation.gnv.artifacts.ressource.RessourceFactory;
 
 import de.intevation.gnv.math.LinearFunction;
 
 import java.text.MessageFormat;
 import java.text.NumberFormat;
+
 import java.util.ArrayList;
 import java.util.Collection;
 import java.util.List;
 import java.util.Locale;
 
+import org.apache.commons.math.FunctionEvaluationException;
+
 import org.apache.commons.math.optimization.OptimizationException;
+
 import org.apache.commons.math.optimization.fitting.CurveFitter;
+
 import org.apache.commons.math.optimization.general.GaussNewtonOptimizer;
-import org.apache.commons.math.FunctionEvaluationException;
 
 import org.apache.log4j.Logger;