Mercurial > dive4elements > gnv-client
view gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java @ 592:3c939c95c477
Removed accidental checkin of Test-Shapefiles
gnv-artifacts/trunk@647 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Tim Englich <tim.englich@intevation.de> |
---|---|
date | Fri, 29 Jan 2010 13:28:28 +0000 |
parents | 44415ae01ddb |
children | 9a828e5a2390 |
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.strtree.STRtree; import java.awt.Dimension; import java.io.Serializable; import java.util.Arrays; import java.util.List; import org.apache.log4j.Logger; /** * @author Sascha L. Teichmann (sascha.teichmann@intevation.de) */ public class Interpolation3D implements Serializable { private static Logger log = Logger.getLogger(Interpolation3D.class); public static final int DEFAULT_WIDTH = 1024; public static final int DEFAULT_HEIGHT = 768; public static final double EPS = 1e-6d; protected int width; protected int height; protected double cellWidth; protected double cellHeight; protected double [] raster; protected double [] depths; public Interpolation3D() { this(DEFAULT_WIDTH, DEFAULT_HEIGHT); } public Interpolation3D(Dimension size) { this(size.width, size.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 getCellWidth() { return cellWidth; } public double getCellHeight() { return cellHeight; } public double [] getRaster() { return raster; } public double [] getDepths() { return depths; } public double getMaxDepth() { double maxDepth = Double.MAX_VALUE; for (int i = depths!=null?depths.length-1:0; i >= 0; --i) { double d = depths[i]; if (!Double.isNaN(d) && d < maxDepth) { maxDepth = d; } } return maxDepth; } public boolean interpolate( List<? extends Coordinate> path, List<? extends XYColumn> 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; } List<GridCell> cells = GridCell.pointsToGridCells( points, Interpolation2D.relevantArea(path, points)); if (cells.isEmpty()) { log.warn("no cells found"); return false; } // put into spatial index to speed up finding neighbors. STRtree spatialIndex = new STRtree(); for (GridCell cell: cells) { spatialIndex.insert(cell.getEnvelope(), cell); } 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; } double cellHeight = Math.abs(maxDepth)/height; if (debug) { log.debug("max depth found: " + maxDepth); log.debug("cell size: " + cellWidth + " x " + cellHeight); } double [] raster = new double[width*height]; Arrays.fill(raster, Double.NaN); Envelope queryBuffer = new Envelope(); GridCell.CellFinder finder = new GridCell.CellFinder(); i = 0; for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) { double depth = depths[i]; if (Double.isNaN(depth) || depth >= 0d) { continue; } linearToMap.locate(p, center); queryBuffer.init( center.x - EPS, center.x + EPS, center.y - EPS, center.y + EPS); finder.prepare(center); spatialIndex.query(queryBuffer, finder); GridCell found = finder.found; if (found == null) { continue; } XYColumn n0 = (XYColumn)found.p1; XYColumn n1 = (XYColumn)found.p2; XYColumn n2 = (XYColumn)found.p3; XYColumn n3 = (XYColumn)found.p4; if (n0.prepare(xyDepth) && n1.prepare(xyDepth) && n2.prepare(xyDepth) && n3.prepare(xyDepth) ) { double y1 = Interpolation2D.interpolate( n0.x, n0.y, n1.x, n1.y, center.x); double y2 = Interpolation2D.interpolate( n2.x, n2.y, n3.x, n3.y, center.x); double z = -cellHeight*0.5; for (int j = i; j < raster.length && z >= depth; z -= cellHeight, j += width) { double v0 = n0.value(z); double v1 = n1.value(z); double v2 = n2.value(z); double v3 = n3.value(z); double z1 = Interpolation2D.interpolate( n0.x, v0, n1.x, v1, center.x); double z2 = Interpolation2D.interpolate( n2.x, v2, n3.x, v3, center.x); double value = Interpolation2D.interpolate( y1, z1, y2, z2, center.y); raster[j] = value; } // XXX: Adjusted depth to create no gap // between last value and ground. depths[i] = z+0.5d*cellHeight; } // down the x/y column } // all along the track this.depths = depths; this.raster = raster; this.cellWidth = cellWidth; this.cellHeight = cellHeight; return true; } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8: