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@474: import com.vividsolutions.jts.index.quadtree.Quadtree;
sascha@474: 
sascha@479: import java.awt.Dimension;
sascha@479: 
sascha@479: import java.io.Serializable;
sascha@479: 
sascha@495: import java.util.ArrayList;
sascha@479: import java.util.Arrays;
sascha@479: import java.util.Collections;
sascha@479: import java.util.List;
sascha@479: 
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<? extends Point2d> points,
sascha@474:         Envelope                boundingBox,
sascha@474:         Dimension               samples,
sascha@474:         XYDepth                 depth
sascha@474:     ) {
sascha@482:         boolean debug = log.isDebugEnabled();
sascha@482: 
sascha@474:         if (points == null || points.isEmpty()) {
sascha@482:             if (debug) {
sascha@482:                 log.debug("no points to interpolate");
sascha@482:             }
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:         Envelope gridEnvelope = new Envelope(boundingBox);
sascha@495: 
sascha@495:         gridEnvelope.expandBy(2d*cellWidth, 2d*cellWidth);
sascha@495: 
sascha@495:         ArrayList<Point2d>relevantPoints =
sascha@495:             new ArrayList<Point2d>();
sascha@495: 
sascha@495:         for (int i = points.size()-1; i >= 0; --i) {
sascha@495:             if (gridEnvelope.contains(points.get(i))) {
sascha@495:                 relevantPoints.add(points.get(i));
sascha@495:             }
sascha@495:         }
sascha@495: 
sascha@501:         double [] buffer = Interpolation2D
sascha@501:             .calculateBufferRemoveOutlier(relevantPoints);
sascha@474:         double dxMax = buffer[0];
sascha@474:         double dyMax = buffer[1];
sascha@474: 
sascha@474:         Quadtree spatialIndex = new Quadtree();
sascha@474: 
sascha@495:         for (int i = relevantPoints.size()-1; i >= 0; --i) {
sascha@495:             Point2d p = relevantPoints.get(i);
sascha@474:             spatialIndex.insert(p.envelope(), p);
sascha@474:         }
sascha@474: 
sascha@482:         if (debug) {
sascha@495:             log.debug("points: "        + relevantPoints.size());
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@495:             log.debug("buffer x: "      + dxMax);
sascha@495:             log.debug("buffer y: "      + dyMax);
sascha@482:         }
sascha@482: 
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@482:         double minX = boundingBox.getMinX();
sascha@482:         double minY = boundingBox.getMinY();
sascha@482: 
sascha@501:         long startTime = System.currentTimeMillis();
sascha@501: 
sascha@498:         int row = 0;
sascha@498:         for (int j = 0; j < H; ++j, row += W) {
sascha@482:             double y = j*cellHeight + 0.5d*cellHeight + minY;
sascha@482:             double x = 0.5d*cellWidth + minX;
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@482:                 if (potential.size() < 4) {
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@501:         long stopTime = System.currentTimeMillis();
sascha@501: 
sascha@501:         if (debug) {
sascha@501:             log.debug("interpolation took: " + (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@474: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :