package de.intevation.gnv.math;

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

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

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

import org.apache.log4j.Logger;

public class Interpolation3D
{
    private static Logger log = Logger.getLogger(Interpolation3D.class);

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

    protected int width;
    protected int height;

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

    public Interpolation3D() {
        this(DEFAULT_WIDTH, DEFAULT_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 [] getRaster() {
        return raster;
    }

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

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

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

        if (debug) {
            log.debug("buffer size: " + dxMax + " / " + dyMax);
        }

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

        for (int i = 0; i < M; ++i) {
            Point2d p = points.get(i);
            spatialIndex.insert(p.envelope(), p);
        }

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

        if (debug) {
            log.debug("max depth found: " + maxDepth);
        }

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

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

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

        Envelope queryBuffer  = new Envelope();
        XYColumn [] neighbors = new XYColumn[4];

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

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

            List potential = spatialIndex.query(queryBuffer);

            L1Comparator invL1 = new L1Comparator(center);
            Collections.sort(potential, invL1);

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

            /* bit code of neighbors
               0---1
               | x |
               2---3
            */

            int mask = 0;
            // reversed order is essential here!
            for (int j = potential.size()-1; j >= 0; --j) {
                XYColumn n = (XYColumn)potential.get(j);
                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
            // and we do not have any gaps.
            if (numNeighbors == 4
            && !neighbors[0].hasIGap(neighbors[1])
            && !neighbors[1].hasJGap(neighbors[3])
            && !neighbors[3].hasIGap(neighbors[2])
            && !neighbors[2].hasJGap(neighbors[0])
            &&  neighbors[0].prepare(xyDepth)
            &&  neighbors[1].prepare(xyDepth)
            &&  neighbors[2].prepare(xyDepth)
            &&  neighbors[3].prepare(xyDepth)
            ) {
                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);
                int j = i;
                for (double z = -cellHeight*0.5; 
                    j < raster.length; z -= cellHeight, j += width) {

                    double v0 = neighbors[0].value(z);
                    double v1 = neighbors[1].value(z);
                    double v2 = neighbors[2].value(z);
                    double v3 = neighbors[3].value(z);

                    double z1 = Interpolation2D.interpolate(
                        neighbors[0].x, v0,
                        neighbors[1].x, v1,
                        center.x);
                    double z2 = Interpolation2D.interpolate(
                        neighbors[2].x, v2,
                        neighbors[3].x, v3,
                        center.x);
                    double value = Interpolation2D.interpolate(
                        y1, z1,
                        y2, z2,
                        center.y);
                    raster[j] = value;
                }
            } // down the x/y column
        } // all along the track

        this.depths = depths;
        this.raster = raster;

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