sascha@424: package de.intevation.gnv.raster;
sascha@424: 
sascha@424: /**
sascha@801:  * A 2D double  valued representation of a raster array.
sascha@780:  * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
sascha@424:  */
sascha@424: public class Raster
sascha@424: {
sascha@801:     /**
sascha@801:      * A kernel can be applied to a raster to fold its values.
sascha@801:      */
sascha@424:     public static class Kernel {
sascha@424: 
sascha@801:         /**
sascha@801:          * The coefficients of the folding matrix.
sascha@801:          */
sascha@424:         public static final class Coeff {
sascha@801:             /**
sascha@801:              * i delta index of this coefficient.
sascha@801:              */
sascha@424:             public int    di;
sascha@801:             /**
sascha@801:              * i delta index of this coefficient.
sascha@801:              */
sascha@424:             public int    dj;
sascha@801:             /**
sascha@801:              * The weight of this coefficient.
sascha@801:              */
sascha@424:             public double c;
sascha@424: 
sascha@801:             /**
sascha@801:              * Default constructor
sascha@801:              */
sascha@424:             public Coeff() {
sascha@424:             }
sascha@424: 
sascha@801:             /**
sascha@801:              * Constructor to set the delta position and the weight of
sascha@801:              * the coefficient.
sascha@801:              * @param di The i delta index.
sascha@801:              * @param dj The j delta index.
sascha@801:              * @param c The weight of the index.
sascha@801:              */
sascha@424:             public Coeff(int di, int dj, double c) {
sascha@424:                 this.di = di;
sascha@424:                 this.dj = dj;
sascha@424:                 this.c  = c;
sascha@424:             }
sascha@424: 
sascha@801:             /**
sascha@801:              * Returns a string representation of this coefficient.
sascha@801:              * @return The string representation.
sascha@801:              */
sascha@801:             @Override
sascha@424:             public String toString() {
sascha@424:                 return String.valueOf(c);
sascha@424:             }
sascha@424: 
sascha@424:         } // class Coeff
sascha@424: 
sascha@801:         /**
sascha@801:          * The coefficients of the folding matrix.
sascha@801:          */
sascha@424:         protected Coeff [] coeffs;
sascha@801:         /**
sascha@801:          * The width of the kernel.
sascha@801:          */
sascha@424:         protected int      width;
sascha@424: 
sascha@801:         /**
sascha@801:          * Default constrcutor.
sascha@801:          */
sascha@424:         public Kernel() {
sascha@424:         }
sascha@424: 
sascha@801:         /**
sascha@801:          * Constructor for a kernel with a given coefficient matrix
sascha@801:          * with a given width.
sascha@801:          * @param coeffs The coefficients.
sascha@801:          * @param width The width of the kernel.
sascha@801:          */
sascha@424:         public Kernel(Coeff [] coeffs, int width) {
sascha@424:             this.coeffs = coeffs;
sascha@424:             this.width  = width;
sascha@424:         }
sascha@424: 
sascha@801:         /**
sascha@801:          * Evaluates the 2D Gauss normal distribution a given x, y
sascha@801:          * and a given sigma.
sascha@801:          * @param x The x location.
sascha@801:          * @param y The y location.
sascha@801:          * @param sigma The sigma of the distribution.
sascha@801:          * @return The value at point x, y.
sascha@801:          */
sascha@424:         public static double gauss(double x, double y, double sigma) {
sascha@424:             double s2 = sigma*sigma;
sascha@424:             return 1.0d/(2.0d*Math.PI*s2)*Math.exp(-(x*x + y*y)/(2.0d*s2));
sascha@424:         }
sascha@424: 
sascha@801:         /**
sascha@801:          * Creates a Gauss kernel with sigma = 1 and width 5.
sascha@801:          * @return The new kernel.
sascha@801:          */
sascha@424:         public static Kernel createGauss() {
sascha@424:             return createGauss(1.0d, 5);
sascha@424:         }
sascha@424: 
sascha@801:         /**
sascha@801:          * Creates a Gauss kernel with a given sigma and width.
sascha@801:          * @param sigma The sigma value
sascha@801:          * @param width  The width of the kernel.
sascha@801:          * @return The new kernel.
sascha@801:          */
sascha@424:         public static Kernel createGauss(double sigma, int width) {
sascha@424:             Coeff [] coeffs = new Coeff[width*width];
sascha@424:             double sum = 0d;
sascha@424:             for (int j = 0; j < width; ++j) {
sascha@424:                 int y = -width/2 + j;
sascha@424:                 for (int i = 0; i < width; ++i) {
sascha@424:                     int x = -width/2 + i;
sascha@424:                     double c = gauss(x, y, sigma);
sascha@424:                     coeffs[j*width + i] = new Coeff(x, y, c);
sascha@424:                     sum += c;
sascha@424:                 }
sascha@424:             }
sascha@424:             sum = 1.0/sum;
sascha@424:             for (int i = 0; i < coeffs.length; ++i) {
sascha@424:                 coeffs[i].c *= sum;
sascha@424:             }
sascha@424:             return new Kernel(coeffs, width);
sascha@424:         }
sascha@424: 
sascha@801:         /**
sascha@801:          * Evaluates this Kernel at the raster at raster index i, j.
sascha@801:          * @param access The accessor to the raster.
sascha@801:          * @param i The i raster index.
sascha@801:          * @param j The j raster index.
sascha@801:          * @return The folded value at index i, j.
sascha@801:          */
sascha@424:         public double fold(Access access, int i, int j) {
sascha@424: 
sascha@424:             double s = 0.0d;
sascha@424: 
sascha@424:             double unused = 0d;
sascha@424: 
sascha@424:             for (int k = 0; k < coeffs.length; ++k) {
sascha@424:                 Coeff coeff = coeffs[k];
sascha@424:                 double v = access.get(i + coeff.di, j + coeff.dj);
sascha@424:                 if (Double.isNaN(v)) {
sascha@424:                     unused += coeff.c;
sascha@424:                 }
sascha@424:                 else {
sascha@424:                     s += v * coeff.c;
sascha@424:                 }
sascha@424:             }
sascha@424: 
sascha@424:             return s*(1.0d - unused);
sascha@424:         }
sascha@424: 
sascha@801:         /**
sascha@801:          * Returns a string representation of this kernel.
sascha@801:          * @return The string representation.
sascha@801:          */
sascha@801:         @Override
sascha@424:         public String toString() {
sascha@424:             StringBuilder sb = new StringBuilder();
sascha@424: 
sascha@424:             int height = coeffs.length/width;
sascha@424: 
sascha@424:             for (int j = 0; j < height; ++j) {
sascha@424:                 if (j > 0) {
sascha@424:                     sb.append(System.getProperty("line.separator"));
sascha@424:                 }
sascha@424:                 for (int i = 0; i < width; ++i) {
sascha@424:                     if (i > 0) {
sascha@424:                         sb.append("\t");
sascha@424:                     }
sascha@424:                     sb.append(coeffs[j*width + i]);
sascha@424:                 }
sascha@424:             }
sascha@424: 
sascha@424:             return sb.toString();
sascha@424:         }
sascha@424:     } // class Kernel
sascha@424: 
sascha@801:     /**
sascha@801:      * Interface of allow special access patterns to the raster,
sascha@801:      * especially at the border.
sascha@801:      */
sascha@424:     public interface Access {
sascha@424: 
sascha@801:         /**
sascha@801:          * Returns the raster value at a given index point.
sascha@801:          * @param i The i index
sascha@801:          * @param j The j index
sascha@801:          * @return The value of the raster.
sascha@801:          */
sascha@424:         double get(int i, int j);
sascha@424: 
sascha@424:     } // interface Access
sascha@424: 
sascha@801:     /**
sascha@801:      * Interface to decouple the lookup which is needed
sascha@801:      * to map double values to integers.
sascha@801:      */
sascha@424:     public interface ValueToIndex {
sascha@424: 
sascha@801:         /**
sascha@801:          * Returns the corresponding integer to a given value.
sascha@801:          * @param value The value to be queried.
sascha@801:          * @return The integer value or -1 if no such value was found.
sascha@801:          */
sascha@424:         int toIndex(double value);
sascha@424: 
sascha@424:     } // interface ValueToIndex;
sascha@424: 
sascha@801:     /**
sascha@801:      * Brings raster to an integer representation with equal sized bins
sascha@801:      * of integer vaules.
sascha@801:      */
sascha@424:     public static class IsoClasses
sascha@424:     implements          ValueToIndex
sascha@424:     {
sascha@801:         /**
sascha@801:          * Linear scaling factor.
sascha@801:          */
sascha@424:         protected double m;
sascha@801:         /**
sascha@801:          * Linear offset.
sascha@801:          */
sascha@424:         protected double b;
sascha@801:         /**
sascha@801:          * Number of classes.
sascha@801:          */
sascha@424:         protected int    numClasses;
sascha@424: 
sascha@801:         /**
sascha@801:          * Constructor to build a lookup to transform
sascha@801:          * an raster into an equals sized bin integer array.
sascha@801:          * @param raster The original raster
sascha@801:          * @param numClasses The number of classes.
sascha@801:          */
sascha@424:         public IsoClasses(Raster raster, int numClasses) {
sascha@424: 
sascha@424:             this.numClasses = numClasses;
sascha@424: 
sascha@424:             double min =  Double.MAX_VALUE;
sascha@424:             double max = -Double.MAX_VALUE;
sascha@424: 
sascha@424:             double [] src = raster.raster;
sascha@424: 
sascha@424:             for (int i = 0; i < src.length; ++i) {
sascha@424:                 double x = src[i];
sascha@424:                 if (!Double.isNaN(x)) {
sascha@424:                     if (x < min) min = x;
sascha@424:                     if (x > max) max = x;
sascha@424:                 }
sascha@424:             }
sascha@424: 
sascha@424:             /* f(min) = 0, f(max) = numClasses - 1
sascha@424: 
sascha@424:                I:  0              = m*min + b
sascha@424:                II: numClasses - 1 = m*max + b
sascha@424: 
sascha@424:                II - I:
sascha@424:                numClasses - 1 = m*(max - min)
sascha@424: 
sascha@424:                m = (numClasses - 1)/(max - min)
sascha@424:                b = m*min
sascha@424:             */
sascha@424: 
sascha@424:             if (max == min) {
sascha@424:                 m = 0;
sascha@424:                 b = 0;
sascha@424:             }
sascha@424:             else {
sascha@424:                 m = (numClasses - 1)/(max - min);
sascha@424:                 b = m*min;
sascha@424:             }
sascha@424:         }
sascha@424: 
sascha@801:         /**
sascha@801:          * Implements the double to integer lookup.
sascha@801:          * @param value The value to map.
sascha@801:          * @return The mapped value.
sascha@801:          */
sascha@424:         public int toIndex(double value) {
sascha@424:             return Double.isNaN(value)
sascha@424:                 ? -1
sascha@424:                 : Math.min(Math.max(0, (int)Math.round(m*value + b)), numClasses-1);
sascha@424:         }
sascha@424:     } // class IsoClasses
sascha@424: 
sascha@801:     /**
sascha@801:      * Internal representation of the raster.
sascha@801:      */
sascha@424:     protected double [] raster;
sascha@801:     /**
sascha@801:      * The width of the raster.
sascha@801:      */
sascha@424:     protected int       width;
sascha@424: 
sascha@801:     /**
sascha@801:      * Default constructor. Creates an 0x0 raster.
sascha@801:      */
sascha@424:     public Raster() {
sascha@424:         this(new double[0], 0);
sascha@424:     }
sascha@424: 
sascha@801:     /**
sascha@801:      * Creates a raster with given values and width.
sascha@801:      * @param raster
sascha@801:      * @param width
sascha@801:      */
sascha@424:     public Raster(double [] raster, int width) {
sascha@424:         this.raster = raster;
sascha@424:         this.width  = width;
sascha@424:     }
sascha@424: 
sascha@801:     /**
sascha@801:      * The width of the raster.
sascha@801:      * @return The width.
sascha@801:      */
sascha@424:     public int getWidth() {
sascha@424:         return width;
sascha@424:     }
sascha@424: 
sascha@801:     /**
sascha@801:      * The height of the raster.
sascha@801:      * @return The height.
sascha@801:      */
sascha@424:     public int getHeight() {
sascha@424:         return raster.length / width;
sascha@424:     }
sascha@424: 
sascha@801:     /**
sascha@801:      * The values of the raster.
sascha@801:      * @return The values.
sascha@801:      */
sascha@495:     public double [] getValues() {
sascha@495:         return raster;
sascha@495:     }
sascha@495: 
sascha@801:     /**
sascha@801:      * Creates an integer representation of this raster using a
sascha@801:      * given value to index transformer.
sascha@801:      * @param valueToIndex The transformer to map double to integer values.
sascha@801:      * @return The new created integer representation.
sascha@801:      */
sascha@424:     public int [] toIndexed(ValueToIndex valueToIndex) {
sascha@424:         int [] dst = new int[raster.length];
sascha@424:         for (int i = 0; i < raster.length; ++i) {
sascha@424:             dst[i] = valueToIndex.toIndex(raster[i]);
sascha@424:         }
sascha@424:         return dst;
sascha@424:     }
sascha@424: 
sascha@801:     /**
sascha@801:      * Creates a new raster by applying a kernel to this raster.
sascha@801:      * Access to the raster is continuous at the border.
sascha@801:      * The original raster is not modified.
sascha@801:      * @param kernel The kernel to be applied.
sascha@801:      * @return The new raster.
sascha@801:      */
sascha@424:     public Raster create(Kernel kernel) {
sascha@424:         double [] dst = new double[raster.length];
sascha@424:         Raster r = new Raster(dst, width);
sascha@424:         r.apply(kernel, continueBorder());
sascha@424:         return r;
sascha@424:     }
sascha@424: 
sascha@801:     /**
sascha@801:      * Creates a new raster by applying a kernel to this raster.
sascha@801:      * Access to the raster is given.
sascha@801:      * @param kernel The kernel to be applied.
sascha@801:      * @param access The access method to the raster.
sascha@801:      */
sascha@424:     public void apply(Kernel kernel, Access access) {
sascha@424:         int height = getHeight();
sascha@424:         for (int j = 0; j < height; ++j) {
sascha@424:             int row = j*width;
sascha@424:             for (int i = 0; i < width; ++i) {
sascha@424:                 raster[row + i] = kernel.fold(access, i, j);
sascha@424:             }
sascha@424:         }
sascha@424:     }
sascha@424: 
sascha@801:     /**
sascha@801:      * Returns a continuous access to the raster. Values
sascha@801:      * at the border are repeated til infinity if indices
sascha@801:      * are accessed outside the defined area.
sascha@801:      * @return The accessor.
sascha@801:      */
sascha@424:     public Access continueBorder() {
sascha@424:         return new Access() {
sascha@424:             int height = getHeight()-1;
sascha@424:             public double get(int i, int j) {
sascha@424:                      if (i < 0)      i = 0;
sascha@424:                 else if (i >= width) i = width-1;
sascha@424:                      if (j < 0)      j = 0;
sascha@424:                 else if (j > height) j = height;
sascha@424:                 return raster[j*width + i];
sascha@424:             }
sascha@424:         };
sascha@424:     }
sascha@424: }
sascha@801: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :