Mercurial > dive4elements > gnv-client
diff gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java @ 474:ab29e4ff2fda
Added area interpolation needed for "Horizontalschnitt"
gnv-artifacts/trunk@540 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Thu, 14 Jan 2010 10:34:05 +0000 |
parents | |
children | d47b478e662b |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gnv-artifacts/src/main/java/de/intevation/gnv/math/AreaInterpolation.java Thu Jan 14 10:34:05 2010 +0000 @@ -0,0 +1,156 @@ +package de.intevation.gnv.math; + +import java.io.Serializable; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.List; +import java.util.HashMap; +import java.util.Collections; + +import java.awt.Dimension; + +import com.vividsolutions.jts.geom.Coordinate; +import com.vividsolutions.jts.geom.Envelope; + +import com.vividsolutions.jts.index.quadtree.Quadtree; + +import org.apache.log4j.Logger; + +/** + * @author Sascha L. Teichmann (sascha.teichmann@intevation.de) + */ +public class AreaInterpolation +implements Serializable +{ + private static Logger log = Logger.getLogger(AreaInterpolation.class); + + protected double [] raster; + + protected int width; + protected int height; + + public AreaInterpolation() { + } + + public int getWidth() { + return width; + } + + public int getHeight() { + return height; + } + + public double [] getRaster() { + return raster; + } + + public boolean interpolate( + List<? extends Point2d> points, + Envelope boundingBox, + Dimension samples, + XYDepth depth + ) { + if (points == null || points.isEmpty()) { + return false; + } + + double [] buffer = Interpolation2D.calculateBuffer(points); + double dxMax = buffer[0]; + double dyMax = buffer[1]; + + Quadtree spatialIndex = new Quadtree(); + + for (int i = points.size()-1; i >= 0; --i) { + Point2d p = points.get(i); + spatialIndex.insert(p.envelope(), p); + } + + int W = samples.width; + int H = samples.height; + + double cellWidth = boundingBox.getWidth() / W; + double cellHeight = boundingBox.getHeight() / H; + + Envelope queryBuffer = new Envelope(); + Point2d [] neighbors = new Point2d[4]; + Coordinate center = new Coordinate(); + L1Comparator invL1 = new L1Comparator(center); + + double [] raster = new double[W*H]; + Arrays.fill(raster, Double.NaN); + + int row = 0; + for (int j = 0; j < H; ++j, row += W) { + double y = j*cellHeight + 0.5d*cellHeight; + double x = 0.5d*cellWidth; + for (int i = 0; i < W; ++i, x += cellWidth) { + queryBuffer.init( + x - dxMax, x + dxMax, + y - dyMax, y + dyMax); + + List potential = spatialIndex.query(queryBuffer); + + if (potential.isEmpty()) { + continue; + } + + center.x = x; center.y = y; + Collections.sort(potential, invL1); + + neighbors[0] = neighbors[1] = + neighbors[2] = neighbors[3] = null; + + int mask = 0; + // reversed order is essential here! + for (int k = potential.size()-1; k >= 0; --k) { + Point2d n = (Point2d)potential.get(k); + int code = n.x > center.x ? 1 : 0; + if (n.y > center.y) code |= 2; + neighbors[code] = n; + mask |= 1 << code; + } + + int numNeighbors = Integer.bitCount(mask); + + // only interpolate if we have all four neighbors + // we do not have any gaps and we are below surface. + if (numNeighbors == 4 + && !neighbors[0].hasIGap(neighbors[1]) + && !neighbors[1].hasJGap(neighbors[3]) + && !neighbors[3].hasIGap(neighbors[2]) + && !neighbors[2].hasJGap(neighbors[0]) + && depth.depth(center) <= 0d + ) { + double z1 = Interpolation2D.interpolate( + neighbors[0].x, neighbors[0].z, + neighbors[1].x, neighbors[1].z, + center.x); + double z2 = Interpolation2D.interpolate( + neighbors[2].x, neighbors[2].z, + neighbors[3].x, neighbors[3].z, + center.x); + double y1 = Interpolation2D.interpolate( + neighbors[0].x, neighbors[0].y, + neighbors[1].x, neighbors[1].y, + center.x); + double y2 = Interpolation2D.interpolate( + neighbors[2].x, neighbors[2].y, + neighbors[3].x, neighbors[3].y, + center.x); + raster[row+i] = Interpolation2D.interpolate( + y1, z1, + y2, z2, + center.y); + } + } + } + + this.raster = raster; + this.width = W; + this.height = H; + + return true; + } +} +// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :