Mercurial > dive4elements > gnv-client
diff geo-backend/src/main/java/de/intevation/gnv/geobackend/sde/datasources/RasterObject.java @ 550:84ba7cbff791
Added bicubic spline interpolation on raster tiles.
geo-backend/trunk@528 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 12 Jan 2010 00:27:53 +0000 |
parents | 0dcf068fb552 |
children | 7615ee5d1345 |
line wrap: on
line diff
--- a/geo-backend/src/main/java/de/intevation/gnv/geobackend/sde/datasources/RasterObject.java Sat Jan 09 16:40:19 2010 +0000 +++ b/geo-backend/src/main/java/de/intevation/gnv/geobackend/sde/datasources/RasterObject.java Tue Jan 12 00:27:53 2010 +0000 @@ -8,6 +8,11 @@ import java.io.Serializable; +import org.apache.commons.math.ArgumentOutsideDomainException; + +import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; + +import org.apache.commons.math.analysis.interpolation.SplineInterpolator; import org.apache.log4j.Logger; @@ -22,6 +27,7 @@ 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; @@ -36,6 +42,9 @@ private double [] rasterData; + private transient PolynomialSplineFunction [] splines; + private transient double [] columns; + public RasterObject() { } @@ -154,6 +163,112 @@ 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; @@ -163,8 +278,9 @@ return Double.NaN; } - if (interpolationType == BILINEAR) { - return interpolateBilinear(px, py); + switch (interpolationType) { + case BILINEAR: return interpolateBilinear(px, py); + case BICUBIC : return interpolateBicubic (px, py); } int posX = Math.min(tileWidth-1, (int)Math.round(px)); @@ -182,6 +298,10 @@ return BILINEAR; } + if ("bicubic".equals(key)) { + return BICUBIC; + } + return NEAREST_NEIGHBOR; } }