ingo@1115: /*
ingo@1115:  * Copyright (c) 2010 by Intevation GmbH
ingo@1115:  *
ingo@1115:  * This program is free software under the LGPL (>=v2.1)
ingo@1115:  * Read the file LGPL.txt coming with the software for details
ingo@1115:  * or visit http://www.gnu.org/licenses/ if it does not exist.
ingo@1115:  */
ingo@1115: 
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@514: import com.vividsolutions.jts.index.strtree.STRtree;
sascha@474: 
sascha@479: import java.awt.Dimension;
sascha@479: 
sascha@479: import java.io.Serializable;
sascha@479: 
sascha@479: import java.util.Arrays;
sascha@479: import java.util.List;
sascha@479: 
sascha@474: import org.apache.log4j.Logger;
sascha@474: 
sascha@474: /**
sascha@808:  * Does an area interpolation for "Horizontalschnitte".
sascha@808:  *
sascha@808:  * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
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@808:     /**
sascha@808:      * The generated raster.
sascha@808:      */
sascha@474:     protected double [] raster;
sascha@474: 
sascha@808:     /**
sascha@808:      * The width of the raster.
sascha@808:      */
sascha@474:     protected int width;
sascha@808:     /**
sascha@808:      * The height of the raster.
sascha@808:      */
sascha@474:     protected int height;
sascha@474: 
sascha@808:     /**
sascha@808:      * Default constructor.
sascha@808:      */
sascha@474:     public AreaInterpolation() {
sascha@474:     }
sascha@474: 
sascha@808:     /**
sascha@808:      * Returns the width of the generated raster.
sascha@808:      * @return the width of the raster.
sascha@808:      */
sascha@474:     public int getWidth() {
sascha@474:         return width;
sascha@474:     }
sascha@474: 
sascha@808:     /**
sascha@808:      * Returns the height of the raster.
sascha@808:      * @return The height of the raster.
sascha@808:      */
sascha@474:     public int getHeight() {
sascha@474:         return height;
sascha@474:     }
sascha@474: 
sascha@808:     /**
sascha@808:      * The generated raster.
ingo@815:      * @return the raster.
sascha@808:      */
sascha@474:     public double [] getRaster() {
sascha@474:         return raster;
sascha@474:     }
sascha@474: 
sascha@808:     /**
sascha@808:      * Epsilon for numerical stability.
sascha@808:      */
sascha@514:     public static final double EPS = 1e-6d;
sascha@514: 
sascha@514: 
sascha@808:     /**
sascha@808:      * Fills a raster by interpolating the input values.
sascha@808:      * @param points The attributed input values.
sascha@808:      * @param boundingBox The world area where to interpolate.
sascha@808:      * @param samples The width and height of the raster.
sascha@808:      * @param depth The callback to query the depth at a given point.
sascha@808:      * @param extrapolationRounds Number of extrapolation point layers
sascha@808:      * to generate before doing the interpolation.
sascha@808:      * @return true if the interpolation succeeds else false.
sascha@808:      */
sascha@474:     public boolean interpolate(
sascha@474:         List<? extends Point2d> points,
sascha@474:         Envelope                boundingBox,
sascha@474:         Dimension               samples,
sascha@593:         XYDepth                 depth,
sascha@593:         int                     extrapolationRounds
sascha@474:     ) {
sascha@482:         boolean debug = log.isDebugEnabled();
sascha@482: 
sascha@474:         if (points == null || points.isEmpty()) {
sascha@514:             log.warn("no points to interpolate");
sascha@514:             return false;
sascha@514:         }
sascha@514: 
sascha@517:         List<GridCell> cells = GridCell.pointsToGridCells(
sascha@528:             points,
sascha@528:             Interpolation2D.relevantArea(
sascha@528:                 boundingBox,
sascha@593:                 points),
sascha@593:             extrapolationRounds);
sascha@514: 
sascha@514:         if (cells.isEmpty()) {
sascha@514:             log.warn("no cells to interpolate");
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: 
sascha@514:         STRtree spatialIndex = new STRtree();
sascha@495: 
sascha@514:         for (GridCell cell: cells) {
sascha@514:             spatialIndex.insert(cell.getEnvelope(), cell);
sascha@474:         }
sascha@474: 
sascha@482:         if (debug) {
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@482:         }
sascha@482: 
sascha@514:         Envelope   queryBuffer     = new Envelope();
sascha@514:         Coordinate center          = new Coordinate();
sascha@514:         GridCell.CellFinder finder = new GridCell.CellFinder();
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@517:         int pos = 0;
sascha@517:         for (int j = 0; j < H; ++j) {
sascha@517: 
sascha@482:             double y = j*cellHeight + 0.5d*cellHeight + minY;
sascha@482:             double x = 0.5d*cellWidth + minX;
sascha@517: 
sascha@517:             for (int end = pos + W; pos < end; ++pos, x += cellWidth) {
sascha@474: 
sascha@514:                 queryBuffer.init(x - EPS, x + EPS, y - EPS, y + EPS);
sascha@514:                 center.x = x; center.y = y;
sascha@514:                 finder.prepare(center);
sascha@514:                 spatialIndex.query(queryBuffer, finder);
sascha@474: 
sascha@514:                 GridCell found = finder.found;
sascha@514: 
sascha@517:                 if (found == null || depth.depth(center) > 0d) {
sascha@474:                     continue;
sascha@474:                 }
sascha@474: 
sascha@514:                 double z1 = Interpolation2D.interpolate(
sascha@514:                     found.p1.x, found.p1.z,
sascha@514:                     found.p2.x, found.p2.z,
sascha@514:                     center.x);
sascha@514:                 double z2 = Interpolation2D.interpolate(
sascha@514:                     found.p3.x, found.p3.z,
sascha@514:                     found.p4.x, found.p4.z,
sascha@514:                     center.x);
sascha@514:                 double y1 = Interpolation2D.interpolate(
sascha@514:                     found.p1.x, found.p1.y,
sascha@514:                     found.p2.x, found.p2.y,
sascha@514:                     center.x);
sascha@514:                 double y2 = Interpolation2D.interpolate(
sascha@514:                     found.p3.x, found.p3.y,
sascha@514:                     found.p4.x, found.p4.y,
sascha@514:                     center.x);
sascha@517:                 raster[pos] = Interpolation2D.interpolate(
sascha@514:                     y1, z1,
sascha@514:                     y2, z2,
sascha@514:                     center.y);
sascha@474:             }
sascha@474:         }
sascha@474: 
sascha@501:         long stopTime = System.currentTimeMillis();
sascha@501: 
sascha@501:         if (debug) {
sascha@801:             log.debug("interpolation took: " +
sascha@801:                 (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@836: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :