/*
 * 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.util.List;

import org.apache.log4j.Logger;

/**
 * Interpolates along a given line string. This used to generate
 * "horizontale Schnittprofile".
 *
 * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
 */
public final class Interpolation2D
{
    private static Logger log = Logger.getLogger(Interpolation2D.class);

    /**
     * If the number of input points are over this threshold
     * some culling strategies are applied to reduce the amount of points.
     */
    public static final int CULL_POINT_THRESHOLD = Integer.getInteger(
        "gnv.interpolation2d.cull.point.threshold", 1000);

    /**
     * Numerical stability epsilon.
     */
    public static final double EPS = 1e-6d;

    /**
     * Callback to send the interpolated values back to the
     * caller without building large temporary strcutures.
     */
    public interface Consumer {
        /**
         * Sends an interpolated point back to the called.
         * The interpolated parameter is stored in the z component.
         * @param point The interpolated point.
         * @param success true if interpolation at the point
         * succeed, else false.
         */
        void interpolated(Coordinate point, boolean success);
    } // interface Consumer

    private Interpolation2D() {
    }

    /**
     * Calculates the relevant area for a given line string.
     * @param path The line string.
     * @param points The sample points.
     * @return The relevant area.
     */
    public static final Envelope relevantArea(
        List<? extends Coordinate> path,
        List<? extends Coordinate> points
    ) {
        return relevantArea(path, points, CULL_POINT_THRESHOLD);
    }

    /**
     * Calculates the relevant area for a given bounding box.
     * @param pathBBox The bounding box.
     * @param points The sample points.
     * @return The relevant area.
     */
    public static final Envelope relevantArea(
        Envelope                   pathBBox,
        List<? extends Coordinate> points
    ) {
        return relevantArea(pathBBox, points, CULL_POINT_THRESHOLD);
    }

    /**
     * Calculates the relevant area for a given bounding box.
     * @param pathBBox The bounding box.
     * @param points The sample points.
     * @param threshold If the number of sample points
     * are below this threshold no relevant area is calculated.
     * @return The relevant area.
     */
    public static final Envelope relevantArea(
        Envelope                   pathBBox,
        List<? extends Coordinate> points,
        int                        threshold
    ) {
        return points.size() < threshold
            ? null
            : relevantArea(
                pathBBox,
                pointsBoundingBox(points));
    }

    /**
     * Calculates the relevant area for a given line string.
     * @param path The line string.
     * @param points The sample points.
     * @param threshold If the number of sample points
     * are below this threshold no relevant area is calculated.
     * @return The relevant area.
     */
    public static final Envelope relevantArea(
        List<? extends Coordinate> path,
        List<? extends Coordinate> points,
        int                        threshold
    ) {
        return points.size() < threshold || path.isEmpty()
            ? null
            : relevantArea(
                pointsBoundingBox(path),
                pointsBoundingBox(points));
    }

    /**
     * Heuristic function to define a 'relevant area'.
     * @param pathBBox The bounding box of the line line string.
     * @param pointsBBox The bounding box of the sample points.
     * @return The relevant area.
     */
    public static final Envelope relevantArea(
        Envelope pathBBox,
        Envelope pointsBBox
    ) {
        double pathArea   = pathBBox.getWidth()*pathBBox.getHeight();
        double pointsArea = pointsBBox.getWidth()*pointsBBox.getHeight();

        if (pathArea > 0.8d*pointsArea) { return null; }

        double nArea = 1.44d * pathArea;
        if (nArea < 0.1d*pointsArea) nArea = 0.1d*pointsArea;
        double w = pathBBox.getWidth();
        double h = pathBBox.getHeight();
        double [] d = solveQuadratic(1d, w+h, pathArea - nArea);

        if (d == null) { return null; }

        double extra = pos(d);

        pathBBox.expandBy(extra);

        return pathBBox;
    }

    /**
     * Solves quadratic equation a*x*x + b*x + c = 0.
     * @param a a-coefficent.
     * @param b b-coefficent
     * @param c c-coefficent.
     * @return the solution of the equation, or null if not solvable.
     */
    public static final double [] solveQuadratic(
        double a, double b, double c
    ) {
        double d = b*b - 4d*a*c;
        if (d < 0d) { return null; }

        d = Math.sqrt(d);
        a = 1d/(2d*a);
        b = -b;

        return new double [] { a*(b + d), a*(b - d) };
    }

    /**
     * Return the element of a two element array which
     * is greater or equal zero.
     * @param x The two values.
     * @return The value which is greater or equal zero.
     */
    public static final double pos(double [] x) {
        return x[0] >= 0 ? x[0] : x[1];
    }


    /**
     * Calculates the bounding box of a given line string.
     * @param path The line string.
     * @return The bounding box.
     */
    public static Envelope pointsBoundingBox(
        List<? extends Coordinate> path
    ) {
        int N = path.size();
        Envelope area = new Envelope(path.get(N-1));

        for (int i = N-2; i >= 0; --i) {
            area.expandToInclude(path.get(i));
        }

        return area;
    }

    /**
     * Interpolates linearly a number of coordinates and parameter values along
     * a given line string path. The results are issued to a consumer.
     * @param path The line string path.
     * @param points The sample points.
     * @param from Start point as a scalar value linear
     * referenced on the line string.
     * @param to End point of  as a scalar value linear
     * referenced on the line string.
     * @param steps Number of points to be interpolated.
     * @param metrics The used metric.
     * @param consumer The callback to retrieve the result points.
     */
    public static void interpolate(
        List<? extends Coordinate> path,
        List<? extends Point2d>    points,
        double                     from,
        double                     to,
        int                        steps,
        Metrics                    metrics,
        Consumer                   consumer
    ) {
        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;
        }

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

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

        // 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 dP = (to - from)/steps;

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

        int missedInterpolations = 0;
        int interpolations       = 0;

        for (double p = from; p <= to; p += dP) {
            if (!linearToMap.locate(p, center)) {
                continue;
            }
            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) {
                consumer.interpolated(center, false);
                ++missedInterpolations;
                continue;
            }

            Point2d n0 = found.p1;
            Point2d n1 = found.p2;
            Point2d n2 = found.p3;
            Point2d n3 = found.p4;

            double z1 = interpolate(
                n0.x, n0.z,
                n1.x, n1.z,
                center.x);
            double z2 = interpolate(
                n2.x, n2.z,
                n3.x, n3.z,
                center.x);
            double y1 = interpolate(
                n0.x, n0.y,
                n1.x, n1.y,
                center.x);
            double y2 = interpolate(
                n2.x, n2.y,
                n3.x, n3.y,
                center.x);
            center.z = interpolate(
                y1, z1,
                y2, z2,
                center.y);
            consumer.interpolated(center, true);
            ++interpolations;
        }

        if (debug) {
            log.debug("interpolations: " +
                interpolations + " / " + missedInterpolations);
        }
    }

    /**
     * Linear interpolate a value between (x1, y1) and (x2, y2) at
     * a given x-value.
     * @param x1 x component of first point.
     * @param y1 y component of first point.
     * @param x2 x component of second point.
     * @param y2 y component of second point.
     * @param x The x value.
     * @return The intepolated result.
     */
    public static final double interpolate(
        double x1, double y1,
        double x2, double y2,
        double x
    ) {
        if (x2 == x1) {
            return (y1 + y2)*0.5d;
        }
        double m = (y2-y1)/(x2-x1);
        double b = y1 - m*x1;
        return m*x + b;
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :