package de.intevation.gnv.raster;

/**
 * @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
 */
public class Raster
{
    public static class Kernel {

        public static final class Coeff {
            public int    di;
            public int    dj;
            public double c;

            public Coeff() {
            }

            public Coeff(int di, int dj, double c) {
                this.di = di;
                this.dj = dj;
                this.c  = c;
            }

            public String toString() {
                return String.valueOf(c);
            }

        } // class Coeff

        protected Coeff [] coeffs;
        protected int      width;

        public Kernel() {
        }

        public Kernel(Coeff [] coeffs, int width) {
            this.coeffs = coeffs;
            this.width  = width;
        }

        public static double gauss(double x, double y, double sigma) {
            double s2 = sigma*sigma;
            return 1.0d/(2.0d*Math.PI*s2)*Math.exp(-(x*x + y*y)/(2.0d*s2));
        }

        public static Kernel createGauss() {
            return createGauss(1.0d, 5);
        }

        public static Kernel createGauss(double sigma, int width) {
            Coeff [] coeffs = new Coeff[width*width];
            double sum = 0d;
            for (int j = 0; j < width; ++j) {
                int y = -width/2 + j;
                for (int i = 0; i < width; ++i) {
                    int x = -width/2 + i;
                    double c = gauss(x, y, sigma);
                    coeffs[j*width + i] = new Coeff(x, y, c);
                    sum += c;
                }
            }
            sum = 1.0/sum;
            for (int i = 0; i < coeffs.length; ++i) {
                coeffs[i].c *= sum;
            }
            return new Kernel(coeffs, width);
        }

        public double fold(Access access, int i, int j) {

            double s = 0.0d;

            double unused = 0d;

            for (int k = 0; k < coeffs.length; ++k) {
                Coeff coeff = coeffs[k];
                double v = access.get(i + coeff.di, j + coeff.dj);
                if (Double.isNaN(v)) {
                    unused += coeff.c;
                }
                else {
                    s += v * coeff.c;
                }
            }

            return s*(1.0d - unused);
        }

        public String toString() {
            StringBuilder sb = new StringBuilder();

            int height = coeffs.length/width;

            for (int j = 0; j < height; ++j) {
                if (j > 0) {
                    sb.append(System.getProperty("line.separator"));
                }
                for (int i = 0; i < width; ++i) {
                    if (i > 0) {
                        sb.append("\t");
                    }
                    sb.append(coeffs[j*width + i]);
                }
            }

            return sb.toString();
        }
    } // class Kernel

    public interface Access {

        double get(int i, int j);

    } // interface Access

    public interface ValueToIndex {

        int toIndex(double value);

    } // interface ValueToIndex;

    public static class IsoClasses
    implements          ValueToIndex
    {
        protected double m;
        protected double b;
        protected int    numClasses;

        public IsoClasses(Raster raster, int numClasses) {

            this.numClasses = numClasses;

            double min =  Double.MAX_VALUE;
            double max = -Double.MAX_VALUE;

            double [] src = raster.raster;

            for (int i = 0; i < src.length; ++i) {
                double x = src[i];
                if (!Double.isNaN(x)) {
                    if (x < min) min = x;
                    if (x > max) max = x;
                }
            }

            /* f(min) = 0, f(max) = numClasses - 1

               I:  0              = m*min + b
               II: numClasses - 1 = m*max + b

               II - I:
               numClasses - 1 = m*(max - min)

               m = (numClasses - 1)/(max - min)
               b = m*min
            */

            if (max == min) {
                m = 0;
                b = 0;
            }
            else {
                m = (numClasses - 1)/(max - min);
                b = m*min;
            }
        }

        public int toIndex(double value) {
            return Double.isNaN(value)
                ? -1
                : Math.min(Math.max(0, (int)Math.round(m*value + b)), numClasses-1);
        }
    } // class IsoClasses

    protected double [] raster;
    protected int       width;

    public Raster() {
        this(new double[0], 0);
    }

    public Raster(double [] raster, int width) {
        this.raster = raster;
        this.width  = width;
    }

    public int getWidth() {
        return width;
    }

    public int getHeight() {
        return raster.length / width;
    }

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

    public int [] toIndexed(ValueToIndex valueToIndex) {
        int [] dst = new int[raster.length];
        for (int i = 0; i < raster.length; ++i) {
            dst[i] = valueToIndex.toIndex(raster[i]);
        }
        return dst;
    }

    public Raster create(Kernel kernel) {
        double [] dst = new double[raster.length];
        Raster r = new Raster(dst, width);
        r.apply(kernel, continueBorder());
        return r;
    }

    public void apply(Kernel kernel, Access access) {
        int height = getHeight();
        for (int j = 0; j < height; ++j) {
            int row = j*width;
            for (int i = 0; i < width; ++i) {
                raster[row + i] = kernel.fold(access, i, j);
            }
        }
    }

    public Access continueBorder() {
        return new Access() {
            int height = getHeight()-1;
            public double get(int i, int j) {
                     if (i < 0)      i = 0;
                else if (i >= width) i = width-1;
                     if (j < 0)      j = 0;
                else if (j > height) j = height;
                return raster[j*width + i];
            }
        };
    }
}
// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8: