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 Sascha L. Teichmann 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 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 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 :