/*
 * 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.geobackend.sde.datasources;

import com.vividsolutions.jts.geom.Coordinate;

import java.io.Serializable;

import org.apache.commons.math.ArgumentOutsideDomainException;

import org.apache.commons.math.analysis.interpolation.SplineInterpolator;

import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;

import org.apache.log4j.Logger;

/**
 * @author <a href="mailto:tim.englich@intevation.de">Tim Englich</a>
 * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
 */
public class RasterObject
implements   Serializable
{
    private static Logger log = Logger.getLogger(RasterObject.class);

    public static final int NEAREST_NEIGHBOR = 0;
    public static final int BILINEAR         = 1;
    public static final int BICUBIC          = 2;

    private double mx;
    private double bx;
    private double my;
    private double by;

    private int columnIndex;
    private int rowIndex;

    private int tileWidth;
    private int tileHeight;

    private double [] rasterData;

    private transient PolynomialSplineFunction [] splines;
    private transient double [] columns;

    public RasterObject() {
    }

    public RasterObject(
        double mx, double bx,
        double my, double by,
        int       columnIndex,
        int       rowIndex,
        double [] rasterData,
        int       tileWidth,
        int       tileHeight
    ){
        this.mx          = mx;
        this.bx          = bx;
        this.my          = my;
        this.by          = by;
        this.columnIndex = columnIndex;
        this.rowIndex    = rowIndex;
        this.rasterData  = rasterData;
        this.tileWidth   = tileWidth;
        this.tileHeight  = tileHeight;
    }

    public boolean contains(Coordinate coordinate) {
        double px = mx*coordinate.x + bx;
        double py = my*coordinate.y + by;
        return px >= 0d && py >= 0d && px < tileWidth && py < tileHeight;
    }

    public final double get(int posX, int posY) {
        return rasterData[posY*tileWidth + posX];
    }

    public int getTileWidth() {
        return tileWidth;
    }

    public int getTileHeight() {
        return tileWidth;
    }

    public int getColumnIndex() {
        return columnIndex;
    }

    public int getRowIndex() {
        return rowIndex;
    }

    public double getValue(Coordinate coordinate) {
        return getValue(coordinate, NEAREST_NEIGHBOR);
    }

    public final double interpolateBilinear(double px, double py) {

        if (px < 0.5d) { // left border
            if (py < 0.5d) { // upper left edge
                return rasterData[0];
            }
            if (py > tileHeight-0.5d) { // lower left edge
                return rasterData[rasterData.length-tileWidth];
            }
            // center left
            int    i = (int)(py -= 0.5d);
            double t = py - i;
            i *= tileWidth;
            return (1d-t)*rasterData[i] + t*rasterData[i+tileWidth];
        }

        if (px > tileWidth-0.5d) { // right border
            if (py < 0.5d) { // upper right edge
                return rasterData[tileWidth-1];
            }
            if (py > tileHeight-0.5d) { // lower right edge
                return rasterData[rasterData.length-1];
            }
            // center right
            int    i = (int)(py -= 0.5d);
            double t = py - i;
            i = i*tileWidth + tileWidth-1;
            return (1d-t)*rasterData[i] + t*rasterData[i+tileWidth];
        }

        if (py < 0.5d) { // top border
            int    i = (int)(px -= 0.5d);
            double t = px - i;
            return (1d-t)*rasterData[i] + t*rasterData[i+1];
        }

        if (py > tileHeight-0.5d) { // bottom border
            int    i = (int)(px -= 0.5d);
            double t = px - i;
            i += rasterData.length-tileWidth;
            return (1d-t)*rasterData[i] + t*rasterData[i+1];
        }

        // center
        int    i  = (int)(py -= 0.5d);
        int    j  = (int)(px -= 0.5d);
        double ti = py - i;
        double tj = px - j;

        int idx = i*tileWidth + j;

        double v1j1 = rasterData[idx];
        double v2j1 = rasterData[idx+1];

        idx += tileWidth;

        double v1j2 = rasterData[idx];
        double v2j2 = rasterData[idx+1];

        double v1 = (1d-tj)*v1j1 + tj*v2j1;
        double v2 = (1d-tj)*v1j2 + tj*v2j2;

        return (1d-ti)*v1 + ti*v2;
    }

    protected PolynomialSplineFunction [] createSplines() {
        PolynomialSplineFunction [] splines =
            new PolynomialSplineFunction[tileHeight];

        double [] x = new double[tileWidth];
        double [] y = new double[tileWidth];

        for (int i = 0; i < tileWidth; ++i) {
            x[i] = i + 0.5d;
        }

        columns = new double[tileHeight];

        SplineInterpolator interpolator = new SplineInterpolator();

        for (int i = 0; i < tileHeight; ++i) {
            int row = i*tileWidth;
            for (int j = 0; j < tileWidth; ++j) {
                y[j] = rasterData[row + j];
            }
            splines[i] = interpolator.interpolate(x, y);
            columns[i] = i + 0.5d;
        }

        return splines;
    }

    protected synchronized PolynomialSplineFunction [] getSplines() {
        if (splines == null) {
            splines = createSplines();
        }
        return splines;
    }

    protected double [] columnValues(double x) {
        PolynomialSplineFunction [] splines = getSplines();
        double [] y = new double[tileHeight];
        try {
            for (int i = 0; i < tileHeight; ++i) {
                y[i] = splines[i].value(x);
            }
        }
        catch (ArgumentOutsideDomainException aode) {
            log.error(aode.getLocalizedMessage(), aode);
        }
        return y;
    }

    private final double evaluateBicubicVertical(double x, double y) {
        try {
            SplineInterpolator interpolator = new SplineInterpolator();
            double [] cs = columnValues(x);
            PolynomialSplineFunction spline = interpolator.interpolate(
                columns, cs);
            return spline.value(y);
        }
        catch (ArgumentOutsideDomainException aode) {
            log.error(aode.getLocalizedMessage(), aode);
        }
        return 0d;
    }

    private final double evaluateBicubicHorizonal(int index, double x) {
        try {
            PolynomialSplineFunction spline = getSplines()[index];
            return spline.value(x);
        }
        catch (ArgumentOutsideDomainException aode) {
            log.error(aode.getLocalizedMessage(), aode);
        }
        return 0d;
    }

    public double interpolateBicubic(double px, double py) {

        if (px <= 0.5d) { // left
            if (py <= 0.5d) { // upper left
                return rasterData[0];
            }
            if (py >= tileHeight-0.5d) {
                return rasterData[rasterData.length-tileWidth];
            }
            return evaluateBicubicVertical(0.5d, py);
        }

        if (px >= tileWidth-0.5d) { // right
            if (py < 0.5d) { // upper right edge
                return rasterData[tileWidth-1];
            }
            if (py > tileHeight-0.5d) { // lower right edge
                return rasterData[rasterData.length-1];
            }
            return evaluateBicubicVertical(tileHeight-0.5d, py);
        }

        if (py <= 0.5d) { // top
            return evaluateBicubicHorizonal(0, px);
        }

        if (py >= tileHeight-0.5d) { // bottom
            return evaluateBicubicHorizonal(tileHeight-1, px);
        }

        return evaluateBicubicVertical(px, py);
    }

    public double getValue(Coordinate coordinate, int interpolationType) {

        double px = mx*coordinate.x + bx;
        double py = my*coordinate.y + by;

        if (px < 0d || py < 0d || px >= tileWidth || py >= tileHeight) {
            return Double.NaN;
        }

        switch (interpolationType) {
            case BILINEAR: return interpolateBilinear(px, py);
            case BICUBIC : return interpolateBicubic (px, py);
        }

        int posX = Math.min(tileWidth-1,  (int)Math.round(px));
        int posY = Math.min(tileHeight-1, (int)Math.round(py));

        return get(posX, posY);
    }

    public static final int getInterpolationType(String key) {
        if (key == null || (key = key.trim().toLowerCase()).length() == 0) {
            return NEAREST_NEIGHBOR;
        }

        if ("bilinear".equals(key)) {
            return BILINEAR;
        }

        if ("bicubic".equals(key)) {
            return BICUBIC;
        }

        return NEAREST_NEIGHBOR;
    }
}