/*
 * Copyright (c) 2010 by Intevation GmbH
 *
 * This program is free software under the LGPL (>=v2.1)
 * Read the file LGPL.txt coming with the software for details
 * or visit http://www.gnu.org/licenses/ if it does not exist.
 */

package de.intevation.gnv.math;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Envelope;

import com.vividsolutions.jts.index.strtree.STRtree;

import java.awt.Dimension;

import java.io.Serializable;

import java.util.Arrays;
import java.util.List;

import org.apache.log4j.Logger;

/**
 * Does an area interpolation for "Horizontalschnitte".
 *
 * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
 */
public class AreaInterpolation
implements   Serializable
{
    private static Logger log = Logger.getLogger(AreaInterpolation.class);

    /**
     * The generated raster.
     */
    protected double [] raster;

    /**
     * The width of the raster.
     */
    protected int width;
    /**
     * The height of the raster.
     */
    protected int height;

    /**
     * Default constructor.
     */
    public AreaInterpolation() {
    }

    /**
     * Returns the width of the generated raster.
     * @return the width of the raster.
     */
    public int getWidth() {
        return width;
    }

    /**
     * Returns the height of the raster.
     * @return The height of the raster.
     */
    public int getHeight() {
        return height;
    }

    /**
     * The generated raster.
     * @return the raster.
     */
    public double [] getRaster() {
        return raster;
    }

    /**
     * Epsilon for numerical stability.
     */
    public static final double EPS = 1e-6d;


    /**
     * Fills a raster by interpolating the input values.
     * @param points The attributed input values.
     * @param boundingBox The world area where to interpolate.
     * @param samples The width and height of the raster.
     * @param depth The callback to query the depth at a given point.
     * @param extrapolationRounds Number of extrapolation point layers
     * to generate before doing the interpolation.
     * @return true if the interpolation succeeds else false.
     */
    public boolean interpolate(
        List<? extends Point2d> points,
        Envelope                boundingBox,
        Dimension               samples,
        XYDepth                 depth,
        int                     extrapolationRounds
    ) {
        boolean debug = log.isDebugEnabled();

        if (points == null || points.isEmpty()) {
            log.warn("no points to interpolate");
            return false;
        }

        List<GridCell> cells = GridCell.pointsToGridCells(
            points,
            Interpolation2D.relevantArea(
                boundingBox,
                points),
            extrapolationRounds);

        if (cells.isEmpty()) {
            log.warn("no cells to interpolate");
            return false;
        }

        int W = samples.width;
        int H = samples.height;

        double cellWidth  = boundingBox.getWidth()  / W;
        double cellHeight = boundingBox.getHeight() / H;

        STRtree spatialIndex = new STRtree();

        for (GridCell cell: cells) {
            spatialIndex.insert(cell.getEnvelope(), cell);
        }

        if (debug) {
            log.debug("width:  "        + boundingBox.getWidth());
            log.debug("height:  "       + boundingBox.getHeight());
            log.debug("sample width:  " + W);
            log.debug("sample height: " + H);
            log.debug("cell width:  "   + cellWidth);
            log.debug("cell height: "   + cellHeight);
        }

        Envelope   queryBuffer     = new Envelope();
        Coordinate center          = new Coordinate();
        GridCell.CellFinder finder = new GridCell.CellFinder();

        double [] raster = new double[W*H];
        Arrays.fill(raster, Double.NaN);

        double minX = boundingBox.getMinX();
        double minY = boundingBox.getMinY();

        long startTime = System.currentTimeMillis();

        int pos = 0;
        for (int j = 0; j < H; ++j) {

            double y = j*cellHeight + 0.5d*cellHeight + minY;
            double x = 0.5d*cellWidth + minX;

            for (int end = pos + W; pos < end; ++pos, x += cellWidth) {

                queryBuffer.init(x - EPS, x + EPS, y - EPS, y + EPS);
                center.x = x; center.y = y;
                finder.prepare(center);
                spatialIndex.query(queryBuffer, finder);

                GridCell found = finder.found;

                if (found == null || depth.depth(center) > 0d) {
                    continue;
                }

                double z1 = Interpolation2D.interpolate(
                    found.p1.x, found.p1.z,
                    found.p2.x, found.p2.z,
                    center.x);
                double z2 = Interpolation2D.interpolate(
                    found.p3.x, found.p3.z,
                    found.p4.x, found.p4.z,
                    center.x);
                double y1 = Interpolation2D.interpolate(
                    found.p1.x, found.p1.y,
                    found.p2.x, found.p2.y,
                    center.x);
                double y2 = Interpolation2D.interpolate(
                    found.p3.x, found.p3.y,
                    found.p4.x, found.p4.y,
                    center.x);
                raster[pos] = Interpolation2D.interpolate(
                    y1, z1,
                    y2, z2,
                    center.y);
            }
        }

        long stopTime = System.currentTimeMillis();

        if (debug) {
            log.debug("interpolation took: " +
                (stopTime - startTime)/1000f + " secs");
        }

        this.raster = raster;
        this.width  = W;
        this.height = H;

        return true;
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :