sascha@434: package de.intevation.gnv.math;
sascha@434:
sascha@434: import com.vividsolutions.jts.geom.Coordinate;
sascha@434: import com.vividsolutions.jts.geom.Envelope;
sascha@434:
sascha@515: import com.vividsolutions.jts.index.strtree.STRtree;
sascha@515:
sascha@515: import java.awt.Dimension;
sascha@515:
sascha@515: import java.io.Serializable;
sascha@515:
sascha@515: import java.util.Arrays;
sascha@515: import java.util.List;
sascha@434:
sascha@434: import org.apache.log4j.Logger;
sascha@434:
sascha@445: /**
sascha@445: * @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
sascha@445: */
sascha@434: public class Interpolation3D
sascha@445: implements Serializable
sascha@434: {
sascha@434: private static Logger log = Logger.getLogger(Interpolation3D.class);
sascha@434:
sascha@434: public static final int DEFAULT_WIDTH = 1024;
sascha@434: public static final int DEFAULT_HEIGHT = 768;
sascha@434:
sascha@515: public static final double EPS = 1e-6d;
sascha@515:
sascha@434: protected int width;
sascha@434: protected int height;
sascha@434:
sascha@521: protected double cellWidth;
sascha@521: protected double cellHeight;
sascha@521:
sascha@434: protected double [] raster;
sascha@434: protected double [] depths;
sascha@434:
sascha@434: public Interpolation3D() {
sascha@434: this(DEFAULT_WIDTH, DEFAULT_HEIGHT);
sascha@434: }
sascha@434:
sascha@445: public Interpolation3D(Dimension size) {
sascha@445: this(size.width, size.height);
sascha@445: }
sascha@445:
sascha@434: public Interpolation3D(int width, int height) {
sascha@434: this.width = width;
sascha@434: this.height = height;
sascha@434: }
sascha@434:
sascha@434: public int getHeight() {
sascha@434: return height;
sascha@434: }
sascha@434:
sascha@434: public int getWidth() {
sascha@434: return width;
sascha@434: }
sascha@434:
sascha@521: public double getCellWidth() {
sascha@521: return cellWidth;
sascha@521: }
sascha@521:
sascha@521: public double getCellHeight() {
sascha@521: return cellHeight;
sascha@521: }
sascha@521:
sascha@434: public double [] getRaster() {
sascha@434: return raster;
sascha@434: }
sascha@434:
sascha@434: public double [] getDepths() {
sascha@434: return depths;
sascha@434: }
sascha@434:
sascha@445: public double getMaxDepth() {
sascha@445: double maxDepth = Double.MAX_VALUE;
sascha@445: for (int i = depths!=null?depths.length-1:0; i >= 0; --i) {
sascha@446: double d = depths[i];
sascha@446: if (!Double.isNaN(d) && d < maxDepth) {
sascha@446: maxDepth = d;
sascha@445: }
sascha@445: }
sascha@445: return maxDepth;
sascha@445: }
sascha@445:
sascha@434: public boolean interpolate(
sascha@434: List<? extends Coordinate> path,
sascha@434: List<? extends XYColumn> points,
sascha@434: double from,
sascha@434: double to,
sascha@434: Metrics metrics,
sascha@434: XYDepth xyDepth
sascha@434: ) {
sascha@434: boolean debug = log.isDebugEnabled();
sascha@434:
sascha@434: int N = path.size();
sascha@434: int M = points.size();
sascha@434:
sascha@434: if (debug) {
sascha@434: log.debug("Size of path: " + N);
sascha@434: log.debug("Size of points: " + M);
sascha@434: }
sascha@434:
sascha@434: if (M < 1 || N < 2) { // nothing to do
sascha@434: return false;
sascha@434: }
sascha@434:
sascha@518: List<GridCell> cells = GridCell.pointsToGridCells(
sascha@528: points, Interpolation2D.relevantArea(path, points));
sascha@434:
sascha@515: if (cells.isEmpty()) {
sascha@515: log.warn("no cells found");
sascha@515: return false;
sascha@434: }
sascha@434:
sascha@434: // put into spatial index to speed up finding neighbors.
sascha@515: STRtree spatialIndex = new STRtree();
sascha@434:
sascha@515: for (GridCell cell: cells) {
sascha@515: spatialIndex.insert(cell.getEnvelope(), cell);
sascha@434: }
sascha@434:
sascha@434: LinearToMap linearToMap = new LinearToMap(
sascha@434: path, from, to, metrics);
sascha@434:
sascha@434: double [] depths = new double[width];
sascha@434: Arrays.fill(depths, Double.NaN);
sascha@434:
sascha@434: double cellWidth = (to - from)/width;
sascha@434:
sascha@434: double maxDepth = Double.MAX_VALUE;
sascha@434:
sascha@434: int i = 0;
sascha@434: Coordinate center = new Coordinate();
sascha@434: for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
sascha@434: if (linearToMap.locate(p, center)) {
sascha@434: double depth = xyDepth.depth(center);
sascha@434: depths[i] = depth;
sascha@434: if (depth < maxDepth) {
sascha@434: maxDepth = depth;
sascha@434: }
sascha@434: }
sascha@434: }
sascha@434:
sascha@434: if (maxDepth == Double.MAX_VALUE) {
sascha@434: log.warn("no depth found -> no interpolation");
sascha@434: return false;
sascha@434: }
sascha@434:
sascha@434: double cellHeight = Math.abs(maxDepth)/height;
sascha@434:
sascha@434: if (debug) {
sascha@445: log.debug("max depth found: " + maxDepth);
sascha@434: log.debug("cell size: " + cellWidth + " x " + cellHeight);
sascha@434: }
sascha@434:
sascha@434: double [] raster = new double[width*height];
sascha@434: Arrays.fill(raster, Double.NaN);
sascha@434:
sascha@515: Envelope queryBuffer = new Envelope();
sascha@515: GridCell.CellFinder finder = new GridCell.CellFinder();
sascha@434:
sascha@434: i = 0;
sascha@434: for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
sascha@434: double depth = depths[i];
sascha@521: if (Double.isNaN(depth) || depth >= 0d) {
sascha@434: continue;
sascha@434: }
sascha@434: linearToMap.locate(p, center);
sascha@434:
sascha@434: queryBuffer.init(
sascha@515: center.x - EPS, center.x + EPS,
sascha@515: center.y - EPS, center.y + EPS);
sascha@434:
sascha@515: finder.prepare(center);
sascha@515: spatialIndex.query(queryBuffer, finder);
sascha@434:
sascha@515: GridCell found = finder.found;
sascha@515:
sascha@515: if (found == null) {
sascha@515: continue;
sascha@434: }
sascha@434:
sascha@515: XYColumn n0 = (XYColumn)found.p1;
sascha@515: XYColumn n1 = (XYColumn)found.p2;
sascha@515: XYColumn n2 = (XYColumn)found.p3;
sascha@515: XYColumn n3 = (XYColumn)found.p4;
sascha@434:
sascha@515: if (n0.prepare(xyDepth)
sascha@515: && n1.prepare(xyDepth)
sascha@515: && n2.prepare(xyDepth)
sascha@515: && n3.prepare(xyDepth)
sascha@434: ) {
sascha@434: double y1 = Interpolation2D.interpolate(
sascha@515: n0.x, n0.y,
sascha@515: n1.x, n1.y,
sascha@434: center.x);
sascha@434: double y2 = Interpolation2D.interpolate(
sascha@515: n2.x, n2.y,
sascha@515: n3.x, n3.y,
sascha@434: center.x);
sascha@521: double z = -cellHeight*0.5;
sascha@521: for (int j = i;
sascha@453: j < raster.length && z >= depth;
sascha@453: z -= cellHeight, j += width) {
sascha@434:
sascha@515: double v0 = n0.value(z);
sascha@515: double v1 = n1.value(z);
sascha@515: double v2 = n2.value(z);
sascha@515: double v3 = n3.value(z);
sascha@434:
sascha@434: double z1 = Interpolation2D.interpolate(
sascha@515: n0.x, v0,
sascha@515: n1.x, v1,
sascha@434: center.x);
sascha@434: double z2 = Interpolation2D.interpolate(
sascha@515: n2.x, v2,
sascha@515: n3.x, v3,
sascha@434: center.x);
sascha@434: double value = Interpolation2D.interpolate(
sascha@434: y1, z1,
sascha@434: y2, z2,
sascha@434: center.y);
sascha@434: raster[j] = value;
sascha@434: }
sascha@521: // XXX: Adjusted depth to create no gap
sascha@521: // between last value and ground.
sascha@521: depths[i] = z+0.5d*cellHeight;
sascha@434: } // down the x/y column
sascha@434: } // all along the track
sascha@434:
sascha@521: this.depths = depths;
sascha@521: this.raster = raster;
sascha@521: this.cellWidth = cellWidth;
sascha@521: this.cellHeight = cellHeight;
sascha@434:
sascha@434: return true;
sascha@434: }
sascha@434: }
sascha@434: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8: