sascha@474: package de.intevation.gnv.math;
sascha@474:
sascha@474: import com.vividsolutions.jts.geom.Coordinate;
sascha@474: import com.vividsolutions.jts.geom.Envelope;
sascha@474:
sascha@514: import com.vividsolutions.jts.index.strtree.STRtree;
sascha@474:
sascha@479: import java.awt.Dimension;
sascha@479:
sascha@479: import java.io.Serializable;
sascha@479:
sascha@479: import java.util.Arrays;
sascha@479: import java.util.List;
sascha@479:
sascha@474: import org.apache.log4j.Logger;
sascha@474:
sascha@474: /**
sascha@808: * Does an area interpolation for "Horizontalschnitte".
sascha@808: *
sascha@808: * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
sascha@474: */
sascha@474: public class AreaInterpolation
sascha@474: implements Serializable
sascha@474: {
sascha@474: private static Logger log = Logger.getLogger(AreaInterpolation.class);
sascha@474:
sascha@808: /**
sascha@808: * The generated raster.
sascha@808: */
sascha@474: protected double [] raster;
sascha@474:
sascha@808: /**
sascha@808: * The width of the raster.
sascha@808: */
sascha@474: protected int width;
sascha@808: /**
sascha@808: * The height of the raster.
sascha@808: */
sascha@474: protected int height;
sascha@474:
sascha@808: /**
sascha@808: * Default constructor.
sascha@808: */
sascha@474: public AreaInterpolation() {
sascha@474: }
sascha@474:
sascha@808: /**
sascha@808: * Returns the width of the generated raster.
sascha@808: * @return the width of the raster.
sascha@808: */
sascha@474: public int getWidth() {
sascha@474: return width;
sascha@474: }
sascha@474:
sascha@808: /**
sascha@808: * Returns the height of the raster.
sascha@808: * @return The height of the raster.
sascha@808: */
sascha@474: public int getHeight() {
sascha@474: return height;
sascha@474: }
sascha@474:
sascha@808: /**
sascha@808: * The generated raster.
ingo@815: * @return the raster.
sascha@808: */
sascha@474: public double [] getRaster() {
sascha@474: return raster;
sascha@474: }
sascha@474:
sascha@808: /**
sascha@808: * Epsilon for numerical stability.
sascha@808: */
sascha@514: public static final double EPS = 1e-6d;
sascha@514:
sascha@514:
sascha@808: /**
sascha@808: * Fills a raster by interpolating the input values.
sascha@808: * @param points The attributed input values.
sascha@808: * @param boundingBox The world area where to interpolate.
sascha@808: * @param samples The width and height of the raster.
sascha@808: * @param depth The callback to query the depth at a given point.
sascha@808: * @param extrapolationRounds Number of extrapolation point layers
sascha@808: * to generate before doing the interpolation.
sascha@808: * @return true if the interpolation succeeds else false.
sascha@808: */
sascha@474: public boolean interpolate(
sascha@474: List<? extends Point2d> points,
sascha@474: Envelope boundingBox,
sascha@474: Dimension samples,
sascha@593: XYDepth depth,
sascha@593: int extrapolationRounds
sascha@474: ) {
sascha@482: boolean debug = log.isDebugEnabled();
sascha@482:
sascha@474: if (points == null || points.isEmpty()) {
sascha@514: log.warn("no points to interpolate");
sascha@514: return false;
sascha@514: }
sascha@514:
sascha@517: List<GridCell> cells = GridCell.pointsToGridCells(
sascha@528: points,
sascha@528: Interpolation2D.relevantArea(
sascha@528: boundingBox,
sascha@593: points),
sascha@593: extrapolationRounds);
sascha@514:
sascha@514: if (cells.isEmpty()) {
sascha@514: log.warn("no cells to interpolate");
sascha@474: return false;
sascha@474: }
sascha@474:
sascha@495: int W = samples.width;
sascha@495: int H = samples.height;
sascha@495:
sascha@495: double cellWidth = boundingBox.getWidth() / W;
sascha@495: double cellHeight = boundingBox.getHeight() / H;
sascha@495:
sascha@514: STRtree spatialIndex = new STRtree();
sascha@495:
sascha@514: for (GridCell cell: cells) {
sascha@514: spatialIndex.insert(cell.getEnvelope(), cell);
sascha@474: }
sascha@474:
sascha@482: if (debug) {
sascha@495: log.debug("width: " + boundingBox.getWidth());
sascha@495: log.debug("height: " + boundingBox.getHeight());
sascha@482: log.debug("sample width: " + W);
sascha@482: log.debug("sample height: " + H);
sascha@495: log.debug("cell width: " + cellWidth);
sascha@495: log.debug("cell height: " + cellHeight);
sascha@482: }
sascha@482:
sascha@514: Envelope queryBuffer = new Envelope();
sascha@514: Coordinate center = new Coordinate();
sascha@514: GridCell.CellFinder finder = new GridCell.CellFinder();
sascha@474:
sascha@474: double [] raster = new double[W*H];
sascha@474: Arrays.fill(raster, Double.NaN);
sascha@474:
sascha@482: double minX = boundingBox.getMinX();
sascha@482: double minY = boundingBox.getMinY();
sascha@482:
sascha@501: long startTime = System.currentTimeMillis();
sascha@501:
sascha@517: int pos = 0;
sascha@517: for (int j = 0; j < H; ++j) {
sascha@517:
sascha@482: double y = j*cellHeight + 0.5d*cellHeight + minY;
sascha@482: double x = 0.5d*cellWidth + minX;
sascha@517:
sascha@517: for (int end = pos + W; pos < end; ++pos, x += cellWidth) {
sascha@474:
sascha@514: queryBuffer.init(x - EPS, x + EPS, y - EPS, y + EPS);
sascha@514: center.x = x; center.y = y;
sascha@514: finder.prepare(center);
sascha@514: spatialIndex.query(queryBuffer, finder);
sascha@474:
sascha@514: GridCell found = finder.found;
sascha@514:
sascha@517: if (found == null || depth.depth(center) > 0d) {
sascha@474: continue;
sascha@474: }
sascha@474:
sascha@514: double z1 = Interpolation2D.interpolate(
sascha@514: found.p1.x, found.p1.z,
sascha@514: found.p2.x, found.p2.z,
sascha@514: center.x);
sascha@514: double z2 = Interpolation2D.interpolate(
sascha@514: found.p3.x, found.p3.z,
sascha@514: found.p4.x, found.p4.z,
sascha@514: center.x);
sascha@514: double y1 = Interpolation2D.interpolate(
sascha@514: found.p1.x, found.p1.y,
sascha@514: found.p2.x, found.p2.y,
sascha@514: center.x);
sascha@514: double y2 = Interpolation2D.interpolate(
sascha@514: found.p3.x, found.p3.y,
sascha@514: found.p4.x, found.p4.y,
sascha@514: center.x);
sascha@517: raster[pos] = Interpolation2D.interpolate(
sascha@514: y1, z1,
sascha@514: y2, z2,
sascha@514: center.y);
sascha@474: }
sascha@474: }
sascha@474:
sascha@501: long stopTime = System.currentTimeMillis();
sascha@501:
sascha@501: if (debug) {
sascha@801: log.debug("interpolation took: " +
sascha@801: (stopTime - startTime)/1000f + " secs");
sascha@501: }
sascha@501:
sascha@474: this.raster = raster;
sascha@474: this.width = W;
sascha@474: this.height = H;
sascha@474:
sascha@474: return true;
sascha@474: }
sascha@474: }
sascha@836: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :