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<? extends Coordinate> path,
sascha@434:         List<? extends XYColumn>   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: