# HG changeset patch # User Sascha L. Teichmann # Date 1261582120 0 # Node ID 0eed5749fd6358b63780a4434c4e4e52ac657606 # Parent 828df3ddb7583451543269689ccf9fdc251edd1c Implemented the raster interpolation for the 'Profilschnitt'. gnv-artifacts/trunk@482 c6561f87-3c4e-4783-a992-168aeb5c3f6f diff -r 828df3ddb758 -r 0eed5749fd63 gnv-artifacts/ChangeLog --- a/gnv-artifacts/ChangeLog Wed Dec 23 12:20:25 2009 +0000 +++ b/gnv-artifacts/ChangeLog Wed Dec 23 15:28:40 2009 +0000 @@ -1,3 +1,29 @@ +2009-12-23 Sascha L. Teichmann + + * src/main/java/de/intevation/gnv/math/Interpolation2D.java: + Use local variable debug instead of asking log.isDebugEnabled() + more than once. + + * src/main/java/de/intevation/gnv/math/XYColumn.java: extrapolate + with boundary values. + + * src/main/java/de/intevation/gnv/math/Interpolation3D.java: + Implements a 3D interpolation called 'Profilschnitt' along a + track similiar to the 'Horizontaler Schnittprofil' which takes + all k layers into account. + + At the interpolated (x, y) points columns of parameter values + from surface to ground are interpolated. To do so the four + next neighbor of that columns are figured out. Four + cubic splines are fitted through these parameter values + of these neighbors. Now its possible to continuous eval + the parameter on each. Every entry in the interpolated column + is interpolated bilinear from the four cubic spline interpolated + neighbor values at the respective depth. + + The result is stored into a double valued raster. NaN values + indicates interpolation gaps. + 2009-12-23 Sascha L. Teichmann * src/main/java/de/intevation/gnv/math/XYColumn.java: Added a method diff -r 828df3ddb758 -r 0eed5749fd63 gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java --- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java Wed Dec 23 12:20:25 2009 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation2D.java Wed Dec 23 15:28:40 2009 +0000 @@ -91,10 +91,12 @@ Metrics metrics, Consumer consumer ) { + boolean debug = log.isDebugEnabled(); + int N = path.size(); int M = points.size(); - if (log.isDebugEnabled()) { + if (debug) { log.debug("Size of path: " + N); log.debug("Size of points: " + M); } @@ -107,7 +109,7 @@ double dxMax = buffer[0]; double dyMax = buffer[1]; - if (log.isDebugEnabled()) { + if (debug) { log.debug("buffer size: " + dxMax + " / " + dyMax); } @@ -204,7 +206,7 @@ } } - if (log.isDebugEnabled()) { + if (debug) { log.debug("interpolations: " + interpolations + " / " + missedInterpolations); } diff -r 828df3ddb758 -r 0eed5749fd63 gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java Wed Dec 23 15:28:40 2009 +0000 @@ -0,0 +1,223 @@ +package de.intevation.gnv.math; + +import java.util.List; +import java.util.Arrays; +import java.util.Collections; + +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.Envelope; + +import com.vividsolutions.jts.index.quadtree.Quadtree; + +import org.apache.log4j.Logger; + +public class Interpolation3D +{ + private static Logger log = Logger.getLogger(Interpolation3D.class); + + public static final int DEFAULT_WIDTH = 1024; + public static final int DEFAULT_HEIGHT = 768; + + protected int width; + protected int height; + + protected double [] raster; + protected double [] depths; + + public Interpolation3D() { + this(DEFAULT_WIDTH, DEFAULT_HEIGHT); + } + + public Interpolation3D(int width, int height) { + this.width = width; + this.height = height; + } + + public int getHeight() { + return height; + } + + public int getWidth() { + return width; + } + + public double [] getRaster() { + return raster; + } + + public double [] getDepths() { + return depths; + } + + public boolean interpolate( + List path, + List points, + double from, + double to, + Metrics metrics, + XYDepth xyDepth + ) { + 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 false; + } + + double [] buffer = Interpolation2D.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 [] depths = new double[width]; + Arrays.fill(depths, Double.NaN); + + double cellWidth = (to - from)/width; + + double maxDepth = Double.MAX_VALUE; + + int i = 0; + Coordinate center = new Coordinate(); + for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) { + if (linearToMap.locate(p, center)) { + double depth = xyDepth.depth(center); + depths[i] = depth; + if (depth < maxDepth) { + maxDepth = depth; + } + } + } + + if (maxDepth == Double.MAX_VALUE) { + log.warn("no depth found -> no interpolation"); + return false; + } + + if (debug) { + log.debug("max depth found: " + maxDepth); + } + + double cellHeight = Math.abs(maxDepth)/height; + + if (debug) { + log.debug("cell size: " + cellWidth + " x " + cellHeight); + } + + double [] raster = new double[width*height]; + Arrays.fill(raster, Double.NaN); + + Envelope queryBuffer = new Envelope(); + XYColumn [] neighbors = new XYColumn[4]; + + i = 0; + for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) { + double depth = depths[i]; + if (Double.isNaN(depth)) { + continue; + } + linearToMap.locate(p, center); + + queryBuffer.init( + center.x - dxMax, center.x + dxMax, + center.y - dyMax, center.y + dyMax); + + List potential = spatialIndex.query(queryBuffer); + + L1Comparator invL1 = new L1Comparator(center); + 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 j = potential.size()-1; j >= 0; --j) { + XYColumn n = (XYColumn)potential.get(j); + 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]) + && neighbors[0].prepare(xyDepth) + && neighbors[1].prepare(xyDepth) + && neighbors[2].prepare(xyDepth) + && neighbors[3].prepare(xyDepth) + ) { + 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); + int j = i; + for (double z = -cellHeight*0.5; + j < raster.length; z -= cellHeight, j += width) { + + double v0 = neighbors[0].value(z); + double v1 = neighbors[1].value(z); + double v2 = neighbors[2].value(z); + double v3 = neighbors[3].value(z); + + double z1 = Interpolation2D.interpolate( + neighbors[0].x, v0, + neighbors[1].x, v1, + center.x); + double z2 = Interpolation2D.interpolate( + neighbors[2].x, v2, + neighbors[3].x, v3, + center.x); + double value = Interpolation2D.interpolate( + y1, z1, + y2, z2, + center.y); + raster[j] = value; + } + } // down the x/y column + } // all along the track + + this.depths = depths; + this.raster = raster; + + return true; + } +} +// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8: diff -r 828df3ddb758 -r 0eed5749fd63 gnv-artifacts/src/main/java/de/intevation/gnv/math/XYColumn.java --- a/gnv-artifacts/src/main/java/de/intevation/gnv/math/XYColumn.java Wed Dec 23 12:20:25 2009 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/XYColumn.java Wed Dec 23 15:28:40 2009 +0000 @@ -24,7 +24,7 @@ extends Point2d implements UnivariateRealFunction { - private static Logger log = Logger.getLogger(Interpolation2D.class); + private static Logger log = Logger.getLogger(XYColumn.class); protected List values; @@ -50,6 +50,11 @@ public double value(double depth) { try { if (curve != null) { + HeightValue h = values.get(0); + // extrapolate beyond boundaries by repeating + if (depth > h.z) return h.v; + h = values.get(values.size()-1); + if (depth < h.z) return h.v; return curve.value(depth); } } @@ -60,7 +65,7 @@ return Double.NaN; } - public void prepare(XYDepth xyDepth) { + public boolean prepare(XYDepth xyDepth) { int N = values.size(); if (curve == null && N > 0) { if (N == 1) { @@ -103,13 +108,16 @@ } catch (MathException me) { log.error("interpolation failed", me); + return false; } } } } else { log.error("no points for interpolation"); + return false; } + return true; } protected UnivariateRealInterpolator getInterpolator() {