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;

/**
 * Interpolates parameter values along a given line string from surface
 * to the sea ground to generate "Profilschnitte".
 *
 * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
 */
public class Interpolation3D
implements   Serializable
{
    private static Logger log = Logger.getLogger(Interpolation3D.class);

    /**
     * The default width of the interpolation: {@value}
     */
    public static final int DEFAULT_WIDTH  = 1024;

    /**
     * The default height of the interpolation: {@value}
     */
    public static final int DEFAULT_HEIGHT =  768;

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

    /**
     * The width of the interpolation.
     */
    protected int width;

    /**
     * The height of the interpolation.
     */
    protected int height;

    /**
     * The cell width of the interpolation in world units.
     */
    protected double cellWidth;

    /**
     * The cell height of the interpolation in world units.
     */
    protected double cellHeight;

    /**
     * The coordinates of the interpolation.
     */
    protected Coordinate[] coordinates;

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

    /**
     * The sea ground depth along the line string.
     */
    protected double [] depths;

    /**
     * Default constructor.
     */
    public Interpolation3D() {
        this(DEFAULT_WIDTH, DEFAULT_HEIGHT);
    }

    /**
     * Constructor to create a Interpolation3D with a given size.
     * @param size The size of the interpolation.
     */
    public Interpolation3D(Dimension size) {
        this(size.width, size.height);
    }

    /**
     * Constructor to create a Interpolation3D with a given size.
     * @param width The width of the interpolation.
     * @param height the height of the interpolation.
     */
    public Interpolation3D(int width, int height) {
        this.width  = width;
        this.height = height;
    }

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

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

    /**
     * Returns the cell width of the interpolation in world units.
     * @return The cell width of the interpolation in world units.
     */
    public double getCellWidth() {
        return cellWidth;
    }

    /**
     * Returns the cell height of the interpolation in world units.
     * @return The cell height of the interpolation in world units.
     */
    public double getCellHeight() {
        return cellHeight;
    }

    /**
     * Returns the coordinates used for the interpolation.
     * @return the coordinates used for the interpolation.
     */
    public Coordinate[] getCoordinates() {
        return coordinates;
    }

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

    /**
     * Returns the depths along the line string path.
     * @return The depth along the line string path.
     */
    public double [] getDepths() {
        return depths;
    }

    /**
     * Returns the deepest depth along the line string path.
     * @return The deepest depth along the line string path.
     */
    public double getMaxDepth() {
        double maxDepth = Double.MAX_VALUE;
        for (int i = depths!=null?depths.length-1:0; i >= 0; --i) {
            double d = depths[i];
            if (!Double.isNaN(d) && d < maxDepth) {
                maxDepth = d;
            }
        }
        return maxDepth;
    }

    /**
     * Interpolates parameters along a given line string path from the surface
     * to the sea ground.
     * @param path The line string path.
     * @param points The sample points.
     * @param from Start point in scalar terms of the line string.
     * @param to End point in scalar terms of the line string.
     * @param metrics The used metric.
     * @param xyDepth The callback to query the depth at a given point.
     * @return true if the interpolation succeeds, else false.
     */
    public boolean interpolate(
        List<? extends Coordinate> path,
        List<? extends XYColumn>   points,
        double                     from,
        double                     to,
        Metrics                    metrics,
        XYDepth                    xyDepth
    ) {
        boolean debug = log.isDebugEnabled();

        int N = path.size();
        int M = points.size();

        if (debug) {
            log.debug("Size of path: " + N);
            log.debug("Size of points: " + M);
        }

        if (M < 1 || N < 2) { // nothing to do
            return false;
        }

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

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

        // put into spatial index to speed up finding neighbors.
        STRtree spatialIndex = new STRtree();

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

        LinearToMap linearToMap = new LinearToMap(
            path, from, to, metrics);

        double [] depths = new double[width];
        Arrays.fill(depths, Double.NaN);

        Coordinate[] coordinates = new Coordinate[width];

        double cellWidth = (to - from)/width;

        double maxDepth = Double.MAX_VALUE;

        int i = 0;
        Coordinate center = new Coordinate();
        for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
            if (linearToMap.locate(p, center)) {
                coordinates[i] = (Coordinate) center.clone();
                double depth = xyDepth.depth(center);
                depths[i] = depth;
                if (depth < maxDepth) {
                    maxDepth = depth;
                }
            }
        }

        if (maxDepth == Double.MAX_VALUE) {
            log.warn("no depth found -> no interpolation");
            return false;
        }

        double cellHeight = Math.abs(maxDepth)/height;

        if (debug) {
            log.debug("max depth found: " + maxDepth);
            log.debug("cell size: " + cellWidth + " x " + cellHeight);
        }

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

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

        i = 0;
        for (double p = cellWidth*0.5; i < depths.length; ++i, p += cellWidth) {
            double depth = depths[i];
            if (Double.isNaN(depth) || depth >= 0d) {
                continue;
            }
            linearToMap.locate(p, center);

            queryBuffer.init(
                center.x - EPS, center.x + EPS,
                center.y - EPS, center.y + EPS);

            finder.prepare(center);
            spatialIndex.query(queryBuffer, finder);

            GridCell found = finder.found;

            if (found == null) {
                continue;
            }

            XYColumn n0 = (XYColumn)found.p1;
            XYColumn n1 = (XYColumn)found.p2;
            XYColumn n2 = (XYColumn)found.p3;
            XYColumn n3 = (XYColumn)found.p4;

            if (n0.prepare(xyDepth)
            &&  n1.prepare(xyDepth)
            &&  n2.prepare(xyDepth)
            &&  n3.prepare(xyDepth)
            ) {
                double y1 = Interpolation2D.interpolate(
                    n0.x, n0.y,
                    n1.x, n1.y,
                    center.x);
                double y2 = Interpolation2D.interpolate(
                    n2.x, n2.y,
                    n3.x, n3.y,
                    center.x);
                double z = -cellHeight*0.5;
                for (int j = i;
                    j < raster.length && z >= depth;
                    z -= cellHeight, j += width) {

                    double v0 = n0.value(z);
                    double v1 = n1.value(z);
                    double v2 = n2.value(z);
                    double v3 = n3.value(z);

                    double z1 = Interpolation2D.interpolate(
                        n0.x, v0,
                        n1.x, v1,
                        center.x);
                    double z2 = Interpolation2D.interpolate(
                        n2.x, v2,
                        n3.x, v3,
                        center.x);
                    double value = Interpolation2D.interpolate(
                        y1, z1,
                        y2, z2,
                        center.y);
                    raster[j] = value;
                }
                // XXX: Adjusted depth to create no gap
                // between last value and ground.
                depths[i] = z+0.5d*cellHeight;
            } // down the x/y column
        } // all along the track

        this.coordinates = coordinates;
        this.depths      = depths;
        this.raster      = raster;
        this.cellWidth   = cellWidth;
        this.cellHeight  = cellHeight;

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