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;

/**
 * @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
 */
public class Interpolation3D
implements   Serializable
{
    private static Logger log = Logger.getLogger(Interpolation3D.class);

    public static final int CULL_POINT_THRESHOLD = Integer.getInteger(
        "gnv.interpolation3d.cull.point.threshold", 1000);

    public static final int DEFAULT_WIDTH  = 1024;
    public static final int DEFAULT_HEIGHT =  768;

    public static final double EPS = 1e-6d;

    protected int width;
    protected int height;

    protected double cellWidth;
    protected double cellHeight;

    protected double [] raster;
    protected double [] depths;

    public Interpolation3D() {
        this(DEFAULT_WIDTH, DEFAULT_HEIGHT);
    }

    public Interpolation3D(Dimension size) {
        this(size.width, size.height);
    }

    public Interpolation3D(int width, int height) {
        this.width  = width;
        this.height = height;
    }

    public int getHeight() {
        return height;
    }

    public int getWidth() {
        return width;
    }

    public double getCellWidth() {
        return cellWidth;
    }

    public double getCellHeight() {
        return cellHeight;
    }

    public double [] getRaster() {
        return raster;
    }

    public double [] getDepths() {
        return depths;
    }

    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;
    }

    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;
        }

        Envelope relevantArea = M > CULL_POINT_THRESHOLD
			? Interpolation2D.pathBoundingBox(path, 0.05d)
			: null;
        
        List<GridCell> cells = GridCell.pointsToGridCells(
            points, relevantArea);

        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);

        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)) {
                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.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: