package de.intevation.gnv.math;

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

import com.vividsolutions.jts.index.quadtree.Quadtree;

import java.awt.Dimension;

import java.io.Serializable;

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

import org.apache.log4j.Logger;

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

    protected double [] raster;

    protected int width;
    protected int height;

    public AreaInterpolation() {
    }

    public int getWidth() {
        return width;
    }

    public int getHeight() {
        return height;
    }

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

    public boolean interpolate(
        List<? extends Point2d> points,
        Envelope                boundingBox,
        Dimension               samples,
        XYDepth                 depth
    ) {
        boolean debug = log.isDebugEnabled();

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

        double [] buffer = Interpolation2D.calculateBuffer(points);
        double dxMax = buffer[0];
        double dyMax = buffer[1];

        Quadtree spatialIndex = new Quadtree();

        for (int i = points.size()-1; i >= 0; --i) {
            Point2d p = points.get(i);
            spatialIndex.insert(p.envelope(), p);
        }

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

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

        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);
            log.debug("buffer x: " + dxMax);
            log.debug("buffer y: " + dyMax);
        }

        Envelope     queryBuffer = new Envelope();
        Point2d []   neighbors   = new Point2d[4];
        Coordinate   center      = new Coordinate();
        L1Comparator invL1       = new L1Comparator(center);

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

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

        int row = 0;
        for (int j = 0; j < H; ++j, row += W) {
            double y = j*cellHeight + 0.5d*cellHeight + minY;
            double x = 0.5d*cellWidth + minX;
            for (int i = 0; i < W; ++i, x += cellWidth) {
                queryBuffer.init(
                    x - dxMax, x + dxMax, 
                    y - dyMax, y + dyMax);

                List potential = spatialIndex.query(queryBuffer);

                if (potential.size() < 4) {
                    continue;
                }

                center.x = x; center.y = y;
                Collections.sort(potential, invL1);

                neighbors[0] = neighbors[1] = 
                neighbors[2] = neighbors[3] = null;

                int mask = 0;
                // reversed order is essential here!
                for (int k = potential.size()-1; k >= 0; --k) {
                    Point2d n = (Point2d)potential.get(k);
                    int code = n.x > center.x ? 1 : 0;
                    if (n.y > center.y) code |= 2;
                    neighbors[code] = n;
                    mask |= 1 << code;
                }

                int numNeighbors = Integer.bitCount(mask);

                // only interpolate if we have all four neighbors
                // we do not have any gaps and we are below surface.
                if (numNeighbors == 4
                && !neighbors[0].hasIGap(neighbors[1])
                && !neighbors[1].hasJGap(neighbors[3])
                && !neighbors[3].hasIGap(neighbors[2])
                && !neighbors[2].hasJGap(neighbors[0])
                && depth.depth(center) <= 0d
                ) {
                    double z1 = Interpolation2D.interpolate(
                        neighbors[0].x, neighbors[0].z,
                        neighbors[1].x, neighbors[1].z,
                        center.x);
                    double z2 = Interpolation2D.interpolate(
                        neighbors[2].x, neighbors[2].z,
                        neighbors[3].x, neighbors[3].z,
                        center.x);
                    double y1 = Interpolation2D.interpolate(
                        neighbors[0].x, neighbors[0].y,
                        neighbors[1].x, neighbors[1].y,
                        center.x);
                    double y2 = Interpolation2D.interpolate(
                        neighbors[2].x, neighbors[2].y,
                        neighbors[3].x, neighbors[3].y,
                        center.x);
                    raster[row+i] = Interpolation2D.interpolate(
                        y1, z1,
                        y2, z2,
                        center.y);
                }
            }
        }

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

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