Mercurial > dive4elements > gnv-client
view gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java @ 505:7ff916744f40
Solved issue152. Time intervals of timeseries axis are defined by JFreeChart automatically.
gnv-artifacts/trunk@588 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Ingo Weinzierl <ingo.weinzierl@intevation.de> |
---|---|
date | Wed, 20 Jan 2010 17:54:48 +0000 |
parents | 70adafe2b9d5 |
children | d9d933e06875 |
line wrap: on
line source
package de.intevation.gnv.math; 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 (sascha.teichmann@intevation.de) */ public final class Interpolation2D { private static Logger log = Logger.getLogger(Interpolation2D.class); public interface Consumer { void interpolated(Coordinate point, boolean success); } // interface Consumer 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 ) { 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); } 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 }; } public static void interpolate( List<? extends Coordinate> path, List<? extends Point2d> points, double from, double to, int steps, Metrics metrics, Consumer consumer ) { boolean debug = log.isDebugEnabled(); int N = path.size(); int M = points.size(); if (debug) { log.debug("Size of path: " + N); log.debug("Size of points: " + M); } if (M < 1 || N < 2) { // nothing to do return; } double [] buffer = calculateBuffer(points); double dxMax = buffer[0]; double dyMax = buffer[1]; if (debug) { log.debug("buffer size: " + dxMax + " / " + dyMax); } // 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]; int missedInterpolations = 0; int interpolations = 0; L1Comparator invL1 = new L1Comparator(center); 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; /* 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[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; } } if (debug) { log.debug("interpolations: " + interpolations + " / " + missedInterpolations); } } 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: