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 Sascha L. Teichmann 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 :