sascha@474: package de.intevation.gnv.math; sascha@474: sascha@474: import java.io.Serializable; sascha@474: sascha@474: import java.util.ArrayList; sascha@474: import java.util.Arrays; sascha@474: import java.util.List; sascha@474: import java.util.HashMap; sascha@474: import java.util.Collections; sascha@474: sascha@474: import java.awt.Dimension; sascha@474: sascha@474: import com.vividsolutions.jts.geom.Coordinate; sascha@474: import com.vividsolutions.jts.geom.Envelope; sascha@474: sascha@474: import com.vividsolutions.jts.index.quadtree.Quadtree; sascha@474: sascha@474: import org.apache.log4j.Logger; sascha@474: sascha@474: /** sascha@474: * @author Sascha L. Teichmann (sascha.teichmann@intevation.de) 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@474: protected double [] raster; sascha@474: sascha@474: protected int width; sascha@474: protected int height; sascha@474: sascha@474: public AreaInterpolation() { sascha@474: } sascha@474: sascha@474: public int getWidth() { sascha@474: return width; sascha@474: } sascha@474: sascha@474: public int getHeight() { sascha@474: return height; sascha@474: } sascha@474: sascha@474: public double [] getRaster() { sascha@474: return raster; sascha@474: } sascha@474: sascha@474: public boolean interpolate( sascha@474: List points, sascha@474: Envelope boundingBox, sascha@474: Dimension samples, sascha@474: XYDepth depth sascha@474: ) { sascha@474: if (points == null || points.isEmpty()) { sascha@474: return false; sascha@474: } sascha@474: sascha@474: double [] buffer = Interpolation2D.calculateBuffer(points); sascha@474: double dxMax = buffer[0]; sascha@474: double dyMax = buffer[1]; sascha@474: sascha@474: Quadtree spatialIndex = new Quadtree(); sascha@474: sascha@474: for (int i = points.size()-1; i >= 0; --i) { sascha@474: Point2d p = points.get(i); sascha@474: spatialIndex.insert(p.envelope(), p); sascha@474: } sascha@474: sascha@474: int W = samples.width; sascha@474: int H = samples.height; sascha@474: sascha@474: double cellWidth = boundingBox.getWidth() / W; sascha@474: double cellHeight = boundingBox.getHeight() / H; sascha@474: sascha@474: Envelope queryBuffer = new Envelope(); sascha@474: Point2d [] neighbors = new Point2d[4]; sascha@474: Coordinate center = new Coordinate(); sascha@474: L1Comparator invL1 = new L1Comparator(center); sascha@474: sascha@474: double [] raster = new double[W*H]; sascha@474: Arrays.fill(raster, Double.NaN); sascha@474: sascha@474: int row = 0; sascha@474: for (int j = 0; j < H; ++j, row += W) { sascha@474: double y = j*cellHeight + 0.5d*cellHeight; sascha@474: double x = 0.5d*cellWidth; sascha@474: for (int i = 0; i < W; ++i, x += cellWidth) { sascha@474: queryBuffer.init( sascha@474: x - dxMax, x + dxMax, sascha@474: y - dyMax, y + dyMax); sascha@474: sascha@474: List potential = spatialIndex.query(queryBuffer); sascha@474: sascha@474: if (potential.isEmpty()) { sascha@474: continue; sascha@474: } sascha@474: sascha@474: center.x = x; center.y = y; sascha@474: Collections.sort(potential, invL1); sascha@474: sascha@474: neighbors[0] = neighbors[1] = sascha@474: neighbors[2] = neighbors[3] = null; sascha@474: sascha@474: int mask = 0; sascha@474: // reversed order is essential here! sascha@474: for (int k = potential.size()-1; k >= 0; --k) { sascha@474: Point2d n = (Point2d)potential.get(k); sascha@474: int code = n.x > center.x ? 1 : 0; sascha@474: if (n.y > center.y) code |= 2; sascha@474: neighbors[code] = n; sascha@474: mask |= 1 << code; sascha@474: } sascha@474: sascha@474: int numNeighbors = Integer.bitCount(mask); sascha@474: sascha@474: // only interpolate if we have all four neighbors sascha@474: // we do not have any gaps and we are below surface. sascha@474: if (numNeighbors == 4 sascha@474: && !neighbors[0].hasIGap(neighbors[1]) sascha@474: && !neighbors[1].hasJGap(neighbors[3]) sascha@474: && !neighbors[3].hasIGap(neighbors[2]) sascha@474: && !neighbors[2].hasJGap(neighbors[0]) sascha@474: && depth.depth(center) <= 0d sascha@474: ) { sascha@474: double z1 = Interpolation2D.interpolate( sascha@474: neighbors[0].x, neighbors[0].z, sascha@474: neighbors[1].x, neighbors[1].z, sascha@474: center.x); sascha@474: double z2 = Interpolation2D.interpolate( sascha@474: neighbors[2].x, neighbors[2].z, sascha@474: neighbors[3].x, neighbors[3].z, sascha@474: center.x); sascha@474: double y1 = Interpolation2D.interpolate( sascha@474: neighbors[0].x, neighbors[0].y, sascha@474: neighbors[1].x, neighbors[1].y, sascha@474: center.x); sascha@474: double y2 = Interpolation2D.interpolate( sascha@474: neighbors[2].x, neighbors[2].y, sascha@474: neighbors[3].x, neighbors[3].y, sascha@474: center.x); sascha@474: raster[row+i] = Interpolation2D.interpolate( sascha@474: y1, z1, sascha@474: y2, z2, sascha@474: center.y); sascha@474: } sascha@474: } sascha@474: } sascha@474: 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@474: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :