Mercurial > dive4elements > gnv-client
view gnv-artifacts/src/main/java/de/intevation/gnv/math/Interpolation3D.java @ 1052:b30361464775
The pdf format depends on the chart's aspect ratio (issue290).
gnv-artifacts/trunk@1126 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Ingo Weinzierl <ingo.weinzierl@intevation.de> |
---|---|
date | Wed, 26 May 2010 06:45:03 +0000 |
parents | 43f3c0cd60f2 |
children | f953c9a559d8 |
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; /** * Interpolates parameter values along a given line string from surface * to the sea ground to generate "Profilschnitte". * * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a> */ public class Interpolation3D implements Serializable { private static Logger log = Logger.getLogger(Interpolation3D.class); /** * The default width of the interpolation: {@value} */ public static final int DEFAULT_WIDTH = 1024; /** * The default height of the interpolation: {@value} */ public static final int DEFAULT_HEIGHT = 768; /** * Epsilon for numerical stability. */ public static final double EPS = 1e-6d; /** * The width of the interpolation. */ protected int width; /** * The height of the interpolation. */ protected int height; /** * The cell width of the interpolation in world units. */ protected double cellWidth; /** * The cell height of the interpolation in world units. */ protected double cellHeight; /** * The coordinates of the interpolation. */ protected Coordinate[] coordinates; /** * The generated raster. */ protected double [] raster; /** * The sea ground depth along the line string. */ protected double [] depths; /** * Default constructor. */ public Interpolation3D() { this(DEFAULT_WIDTH, DEFAULT_HEIGHT); } /** * Constructor to create a Interpolation3D with a given size. * @param size The size of the interpolation. */ public Interpolation3D(Dimension size) { this(size.width, size.height); } /** * Constructor to create a Interpolation3D with a given size. * @param width The width of the interpolation. * @param height the height of the interpolation. */ public Interpolation3D(int width, int height) { this.width = width; this.height = height; } /** * Returns the raster height of the interpolation. * @return The raster height of the interpolation. */ public int getHeight() { return height; } /** * Returns the raster width of the interpolation. * @return The raster width of the interpolation. */ public int getWidth() { return width; } /** * Returns the cell width of the interpolation in world units. * @return The cell width of the interpolation in world units. */ public double getCellWidth() { return cellWidth; } /** * Returns the cell height of the interpolation in world units. * @return The cell height of the interpolation in world units. */ public double getCellHeight() { return cellHeight; } /** * Returns the coordinates used for the interpolation. * @return the coordinates used for the interpolation. */ public Coordinate[] getCoordinates() { return coordinates; } /** * Returns the generated raster. * @return The generated raster. */ public double [] getRaster() { return raster; } /** * Returns the depths along the line string path. * @return The depth along the line string path. */ public double [] getDepths() { return depths; } /** * Returns the deepest depth along the line string path. * @return The deepest depth along the line string path. */ 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; } /** * Interpolates parameters along a given line string path from the surface * to the sea ground. * @param path The line string path. * @param points The sample points. * @param from Start point in scalar terms of the line string. * @param to End point in scalar terms of the line string. * @param metrics The used metric. * @param xyDepth The callback to query the depth at a given point. * @return true if the interpolation succeeds, else false. */ 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); Coordinate[] coordinates = new Coordinate[width]; 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)) { coordinates[i] = (Coordinate) center.clone(); 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.coordinates = coordinates; 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 :