Mercurial > dive4elements > gnv-client
changeset 519:4e347624ee7c
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
| author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
|---|---|
| date | Sat, 23 Jan 2010 21:16:45 +0000 |
| parents | 464e03bf786b |
| children | a8f6ca59b26e |
| files | gnv-artifacts/ChangeLog 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/Interpolation3D.java gnv-artifacts/src/main/java/de/intevation/gnv/state/profile/horizontal/HorizontalProfileMeshCrossOutputState.java gnv-artifacts/src/main/java/de/intevation/gnv/utils/WKTUtils.java |
| diffstat | 6 files changed, 137 insertions(+), 277 deletions(-) [+] |
line wrap: on
line diff
--- 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:
--- 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; } }
--- 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) {
--- 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);
--- 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,
--- 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;