sascha@424: package de.intevation.gnv.raster;
sascha@424:
sascha@424: /**
sascha@424: * @author Sascha L. Teichmann (sascha.teichmann@intevation.de)
sascha@424: */
sascha@424: public class Raster
sascha@424: {
sascha@424: public static class Kernel {
sascha@424:
sascha@424: public static final class Coeff {
sascha@424: public int di;
sascha@424: public int dj;
sascha@424: public double c;
sascha@424:
sascha@424: public Coeff() {
sascha@424: }
sascha@424:
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@424: public String toString() {
sascha@424: return String.valueOf(c);
sascha@424: }
sascha@424:
sascha@424: } // class Coeff
sascha@424:
sascha@424: protected Coeff [] coeffs;
sascha@424: protected int width;
sascha@424:
sascha@424: public Kernel() {
sascha@424: }
sascha@424:
sascha@424: public Kernel(Coeff [] coeffs, int width) {
sascha@424: this.coeffs = coeffs;
sascha@424: this.width = width;
sascha@424: }
sascha@424:
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@424: public static Kernel createGauss() {
sascha@424: return createGauss(1.0d, 5);
sascha@424: }
sascha@424:
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@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@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@424: public interface Access {
sascha@424:
sascha@424: double get(int i, int j);
sascha@424:
sascha@424: } // interface Access
sascha@424:
sascha@424: public interface ValueToIndex {
sascha@424:
sascha@424: int toIndex(double value);
sascha@424:
sascha@424: } // interface ValueToIndex;
sascha@424:
sascha@424: public static class IsoClasses
sascha@424: implements ValueToIndex
sascha@424: {
sascha@424: protected double m;
sascha@424: protected double b;
sascha@424: protected int numClasses;
sascha@424:
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@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@424: protected double [] raster;
sascha@424: protected int width;
sascha@424:
sascha@424: public Raster() {
sascha@424: this(new double[0], 0);
sascha@424: }
sascha@424:
sascha@424: public Raster(double [] raster, int width) {
sascha@424: this.raster = raster;
sascha@424: this.width = width;
sascha@424: }
sascha@424:
sascha@424: public int getWidth() {
sascha@424: return width;
sascha@424: }
sascha@424:
sascha@424: public int getHeight() {
sascha@424: return raster.length / width;
sascha@424: }
sascha@424:
sascha@495: public double [] getValues() {
sascha@495: return raster;
sascha@495: }
sascha@495:
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@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@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@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@424: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8: