# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1263998850 0
# Node ID 70adafe2b9d575833d848822ff7c1f76e0e21c4a
# Parent  ca5162aa644d6b072cdff4490bc446d6f342809d
Speed up the "Horizontalschnitte"


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

diff -r ca5162aa644d -r 70adafe2b9d5 gnv-artifacts/ChangeLog
--- a/gnv-artifacts/ChangeLog	Wed Jan 20 08:50:10 2010 +0000
+++ b/gnv-artifacts/ChangeLog	Wed Jan 20 14:47:30 2010 +0000
@@ -1,3 +1,23 @@
+2010-01-20	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
+
+	* doc/conf/conf.xml: Set number of per axis samples to 1024
+	  because generation of "Horizontalschnitte" is much faster
+	  now (60x).
+
+	* src/main/java/de/intevation/gnv/math/Interpolation2D.java:
+	  Added some kind of outlier test when guessing the buffer size
+	  of the spatial index. The speed problem arose from the fact
+	  that to much points are assumed to be neighbors of a given
+	  point. Long distances which differ more than 40% from the 
+	  standard derivation are assumed to be outliers.
+
+	* src/main/java/de/intevation/gnv/math/AreaInterpolation.java: Uses
+	  the outlier aware buffer size guessing now.
+
+	* src/main/java/de/intevation/gnv/raster/JTSMultiPolygonProducer.java,
+	  src/main/java/de/intevation/gnv/raster/JTSMultiLineStringProducer.java:
+	  Removed needless imports.
+
 2010-01-20  Ingo Weinzierl <ingo.weinzierl@intevation.de>
 
 	  Issue148
diff -r ca5162aa644d -r 70adafe2b9d5 gnv-artifacts/doc/conf/conf.xml
--- a/gnv-artifacts/doc/conf/conf.xml	Wed Jan 20 08:50:10 2010 +0000
+++ b/gnv-artifacts/doc/conf/conf.xml	Wed Jan 20 14:47:30 2010 +0000
@@ -420,7 +420,7 @@
 
         <horizontal-cross-section>
             <!-- This section configures the HorizontalCrossSection ("Horizontalschnitt") -->
-            <samples number="256"/>
+            <samples number="1024"/>
             <ground interpolation="bilinear" />
             <result-shapefile-directory path="${artifacts.config.dir}/../shapefiles/"/>
         </horizontal-cross-section>
diff -r ca5162aa644d -r 70adafe2b9d5 gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java	Wed Jan 20 08:50:10 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java	Wed Jan 20 14:47:30 2010 +0000
@@ -77,7 +77,8 @@
             }
         }
 
-        double [] buffer = Interpolation2D.calculateBuffer(relevantPoints);
+        double [] buffer = Interpolation2D
+            .calculateBufferRemoveOutlier(relevantPoints);
         double dxMax = buffer[0];
         double dyMax = buffer[1];
 
@@ -111,6 +112,8 @@
         double minX = boundingBox.getMinX();
         double minY = boundingBox.getMinY();
 
+        long startTime = System.currentTimeMillis();
+
         int row = 0;
         for (int j = 0; j < H; ++j, row += W) {
             double y = j*cellHeight + 0.5d*cellHeight + minY;
@@ -177,6 +180,12 @@
             }
         }
 
+        long stopTime = System.currentTimeMillis();
+
+        if (debug) {
+            log.debug("interpolation took: " + (stopTime - startTime)/1000f + " secs");
+        }
+
         this.raster = raster;
         this.width  = W;
         this.height = H;
diff -r ca5162aa644d -r 70adafe2b9d5 gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Wed Jan 20 08:50:10 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java	Wed Jan 20 14:47:30 2010 +0000
@@ -1,19 +1,24 @@
 package de.intevation.gnv.math;
 
-import java.util.ArrayList;
-import java.util.List;
-import java.util.HashMap;
-import java.util.Collections;
-
 import com.vividsolutions.jts.geom.Coordinate;
 import com.vividsolutions.jts.geom.Envelope;
 
 import com.vividsolutions.jts.index.quadtree.Quadtree;
 
+import gnu.trove.TDoubleArrayList;
+
+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;
 
 /**
- *  @author Sascha L. Teichmann
+ *  @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
  */
 public final class Interpolation2D
 {
@@ -26,6 +31,108 @@
     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
+    ) {
+        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);
+
+            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);
+                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);
+        }
+
+        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);
+                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
     ) {
diff -r ca5162aa644d -r 70adafe2b9d5 gnv-artifacts/src/main/java/de/intevation/gnv/raster/JTSMultiLineStringProducer.java
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/raster/JTSMultiLineStringProducer.java	Wed Jan 20 08:50:10 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/raster/JTSMultiLineStringProducer.java	Wed Jan 20 14:47:30 2010 +0000
@@ -2,11 +2,10 @@
 
 import com.vividsolutions.jts.geom.Coordinate;
 import com.vividsolutions.jts.geom.Geometry;
-import com.vividsolutions.jts.geom.GeometryCollection;
 import com.vividsolutions.jts.geom.GeometryFactory;
-import com.vividsolutions.jts.geom.Polygon;
 import com.vividsolutions.jts.geom.LineString;
 import com.vividsolutions.jts.geom.MultiLineString;
+import com.vividsolutions.jts.geom.Polygon;
 import com.vividsolutions.jts.geom.TopologyException;
 
 import de.intevation.gnv.math.IJKey;
diff -r ca5162aa644d -r 70adafe2b9d5 gnv-artifacts/src/main/java/de/intevation/gnv/raster/JTSMultiPolygonProducer.java
--- a/gnv-artifacts/src/main/java/de/intevation/gnv/raster/JTSMultiPolygonProducer.java	Wed Jan 20 08:50:10 2010 +0000
+++ b/gnv-artifacts/src/main/java/de/intevation/gnv/raster/JTSMultiPolygonProducer.java	Wed Jan 20 14:47:30 2010 +0000
@@ -4,7 +4,6 @@
 
 import com.vividsolutions.jts.geom.Coordinate;
 import com.vividsolutions.jts.geom.Geometry;
-import com.vividsolutions.jts.geom.GeometryCollection;
 import com.vividsolutions.jts.geom.GeometryFactory;
 import com.vividsolutions.jts.geom.LinearRing;
 import com.vividsolutions.jts.geom.MultiPolygon;