sascha@434: package de.intevation.gnv.math; sascha@434: sascha@434: import java.util.List; sascha@434: import java.util.Arrays; sascha@434: import java.util.Collections; sascha@434: sascha@445: import java.io.Serializable; sascha@445: sascha@445: import java.awt.Dimension; sascha@445: sascha@434: import com.vividsolutions.jts.geom.Coordinate; sascha@434: import com.vividsolutions.jts.geom.Envelope; sascha@434: sascha@434: import com.vividsolutions.jts.index.quadtree.Quadtree; 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@434: protected int width; sascha@434: protected int height; sascha@434: 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@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 path, sascha@434: List 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@434: double [] buffer = Interpolation2D.calculateBuffer(points); sascha@434: double dxMax = buffer[0]; sascha@434: double dyMax = buffer[1]; sascha@434: sascha@434: if (debug) { sascha@434: log.debug("buffer size: " + dxMax + " / " + dyMax); sascha@434: } sascha@434: sascha@434: // put into spatial index to speed up finding neighbors. sascha@434: Quadtree spatialIndex = new Quadtree(); sascha@434: sascha@434: for (int i = 0; i < M; ++i) { sascha@434: Point2d p = points.get(i); sascha@434: spatialIndex.insert(p.envelope(), p); 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@434: Envelope queryBuffer = new Envelope(); sascha@434: XYColumn [] neighbors = new XYColumn[4]; 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@434: if (Double.isNaN(depth)) { sascha@434: continue; sascha@434: } sascha@434: linearToMap.locate(p, center); sascha@434: sascha@434: queryBuffer.init( sascha@434: center.x - dxMax, center.x + dxMax, sascha@434: center.y - dyMax, center.y + dyMax); sascha@434: sascha@434: List potential = spatialIndex.query(queryBuffer); sascha@434: sascha@434: L1Comparator invL1 = new L1Comparator(center); sascha@434: Collections.sort(potential, invL1); sascha@434: sascha@434: neighbors[0] = neighbors[1] = sascha@434: neighbors[2] = neighbors[3] = null; sascha@434: sascha@434: /* bit code of neighbors sascha@434: 0---1 sascha@434: | x | sascha@434: 2---3 sascha@434: */ sascha@434: sascha@434: int mask = 0; sascha@434: // reversed order is essential here! sascha@434: for (int j = potential.size()-1; j >= 0; --j) { sascha@434: XYColumn n = (XYColumn)potential.get(j); sascha@434: int code = n.x > center.x ? 1 : 0; sascha@434: if (n.y > center.y) code |= 2; sascha@434: neighbors[code] = n; sascha@434: mask |= 1 << code; sascha@434: } sascha@434: sascha@434: int numNeighbors = Integer.bitCount(mask); sascha@434: sascha@434: // only interpolate if we have all four neighbors sascha@434: // and we do not have any gaps. sascha@434: if (numNeighbors == 4 sascha@434: && !neighbors[0].hasIGap(neighbors[1]) sascha@434: && !neighbors[1].hasJGap(neighbors[3]) sascha@434: && !neighbors[3].hasIGap(neighbors[2]) sascha@434: && !neighbors[2].hasJGap(neighbors[0]) sascha@434: && neighbors[0].prepare(xyDepth) sascha@434: && neighbors[1].prepare(xyDepth) sascha@434: && neighbors[2].prepare(xyDepth) sascha@434: && neighbors[3].prepare(xyDepth) sascha@434: ) { sascha@434: double y1 = Interpolation2D.interpolate( sascha@434: neighbors[0].x, neighbors[0].y, sascha@434: neighbors[1].x, neighbors[1].y, sascha@434: center.x); sascha@434: double y2 = Interpolation2D.interpolate( sascha@434: neighbors[2].x, neighbors[2].y, sascha@434: neighbors[3].x, neighbors[3].y, sascha@434: center.x); sascha@434: int j = i; sascha@434: for (double z = -cellHeight*0.5; sascha@453: j < raster.length && z >= depth; sascha@453: z -= cellHeight, j += width) { sascha@434: sascha@434: double v0 = neighbors[0].value(z); sascha@434: double v1 = neighbors[1].value(z); sascha@434: double v2 = neighbors[2].value(z); sascha@434: double v3 = neighbors[3].value(z); sascha@434: sascha@434: double z1 = Interpolation2D.interpolate( sascha@434: neighbors[0].x, v0, sascha@434: neighbors[1].x, v1, sascha@434: center.x); sascha@434: double z2 = Interpolation2D.interpolate( sascha@434: neighbors[2].x, v2, sascha@434: neighbors[3].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@434: } // down the x/y column sascha@434: } // all along the track sascha@434: sascha@434: this.depths = depths; sascha@434: this.raster = raster; 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: