ingo@1127: /*
ingo@1127:  * Copyright (c) 2010 by Intevation GmbH
ingo@1127:  *
ingo@1127:  * This program is free software under the LGPL (>=v2.1)
ingo@1127:  * Read the file LGPL.txt coming with the software for details
ingo@1127:  * or visit http://www.gnu.org/licenses/ if it does not exist.
ingo@1127:  */
ingo@1127: 
tim@544: package de.intevation.gnv.geobackend.sde.datasources;
tim@544: 
tim@544: import com.vividsolutions.jts.geom.Coordinate;
tim@544: 
sascha@545: import java.io.Serializable;
sascha@545: 
sascha@550: import org.apache.commons.math.ArgumentOutsideDomainException;
sascha@550: 
sascha@552: import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
sascha@552: 
sascha@550: import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
sascha@550: 
sascha@547: import org.apache.log4j.Logger;
sascha@547: 
tim@544: /**
sascha@887:  * @author <a href="mailto:tim.englich@intevation.de">Tim Englich</a>
sascha@887:  * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a>
tim@544:  */
sascha@885: public class RasterObject
sascha@545: implements   Serializable
sascha@545: {
sascha@547:     private static Logger log = Logger.getLogger(RasterObject.class);
sascha@547: 
sascha@547:     public static final int NEAREST_NEIGHBOR = 0;
sascha@547:     public static final int BILINEAR         = 1;
sascha@550:     public static final int BICUBIC          = 2;
sascha@547: 
sascha@547:     private double mx;
sascha@547:     private double bx;
sascha@547:     private double my;
sascha@547:     private double by;
sascha@547: 
tim@544:     private int columnIndex;
sascha@547:     private int rowIndex;
sascha@547: 
sascha@547:     private int tileWidth;
sascha@547:     private int tileHeight;
sascha@547: 
sascha@547:     private double [] rasterData;
sascha@547: 
sascha@550:     private transient PolynomialSplineFunction [] splines;
sascha@550:     private transient double [] columns;
sascha@550: 
sascha@547:     public RasterObject() {
sascha@547:     }
sascha@547: 
sascha@547:     public RasterObject(
sascha@547:         double mx, double bx,
sascha@547:         double my, double by,
sascha@885:         int       columnIndex,
sascha@547:         int       rowIndex,
sascha@547:         double [] rasterData,
sascha@547:         int       tileWidth,
sascha@547:         int       tileHeight
sascha@547:     ){
sascha@547:         this.mx          = mx;
sascha@547:         this.bx          = bx;
sascha@547:         this.my          = my;
sascha@547:         this.by          = by;
tim@544:         this.columnIndex = columnIndex;
sascha@547:         this.rowIndex    = rowIndex;
sascha@547:         this.rasterData  = rasterData;
sascha@547:         this.tileWidth   = tileWidth;
sascha@547:         this.tileHeight  = tileHeight;
tim@544:     }
sascha@547: 
sascha@547:     public boolean contains(Coordinate coordinate) {
sascha@547:         double px = mx*coordinate.x + bx;
sascha@547:         double py = my*coordinate.y + by;
sascha@547:         return px >= 0d && py >= 0d && px < tileWidth && py < tileHeight;
sascha@547:     }
sascha@547: 
sascha@547:     public final double get(int posX, int posY) {
sascha@547:         return rasterData[posY*tileWidth + posX];
sascha@547:     }
sascha@547: 
sascha@547:     public int getTileWidth() {
sascha@547:         return tileWidth;
sascha@547:     }
sascha@547: 
sascha@547:     public int getTileHeight() {
sascha@547:         return tileWidth;
sascha@547:     }
sascha@547: 
sascha@547:     public int getColumnIndex() {
sascha@547:         return columnIndex;
sascha@547:     }
sascha@547: 
sascha@547:     public int getRowIndex() {
sascha@547:         return rowIndex;
tim@544:     }
sascha@885: 
sascha@547:     public double getValue(Coordinate coordinate) {
sascha@547:         return getValue(coordinate, NEAREST_NEIGHBOR);
tim@544:     }
tim@544: 
sascha@548:     public final double interpolateBilinear(double px, double py) {
sascha@548: 
sascha@548:         if (px < 0.5d) { // left border
sascha@548:             if (py < 0.5d) { // upper left edge
sascha@548:                 return rasterData[0];
sascha@548:             }
sascha@548:             if (py > tileHeight-0.5d) { // lower left edge
sascha@548:                 return rasterData[rasterData.length-tileWidth];
sascha@548:             }
sascha@548:             // center left
sascha@548:             int    i = (int)(py -= 0.5d);
sascha@548:             double t = py - i;
sascha@548:             i *= tileWidth;
sascha@548:             return (1d-t)*rasterData[i] + t*rasterData[i+tileWidth];
sascha@548:         }
sascha@548: 
sascha@548:         if (px > tileWidth-0.5d) { // right border
sascha@548:             if (py < 0.5d) { // upper right edge
sascha@548:                 return rasterData[tileWidth-1];
sascha@548:             }
sascha@548:             if (py > tileHeight-0.5d) { // lower right edge
sascha@548:                 return rasterData[rasterData.length-1];
sascha@548:             }
sascha@548:             // center right
sascha@548:             int    i = (int)(py -= 0.5d);
sascha@548:             double t = py - i;
sascha@548:             i = i*tileWidth + tileWidth-1;
sascha@548:             return (1d-t)*rasterData[i] + t*rasterData[i+tileWidth];
sascha@548:         }
sascha@548: 
sascha@548:         if (py < 0.5d) { // top border
sascha@548:             int    i = (int)(px -= 0.5d);
sascha@548:             double t = px - i;
sascha@548:             return (1d-t)*rasterData[i] + t*rasterData[i+1];
sascha@548:         }
sascha@548: 
sascha@548:         if (py > tileHeight-0.5d) { // bottom border
sascha@548:             int    i = (int)(px -= 0.5d);
sascha@548:             double t = px - i;
sascha@548:             i += rasterData.length-tileWidth;
sascha@548:             return (1d-t)*rasterData[i] + t*rasterData[i+1];
sascha@548:         }
sascha@548: 
sascha@548:         // center
sascha@548:         int    i  = (int)(py -= 0.5d);
sascha@548:         int    j  = (int)(px -= 0.5d);
sascha@548:         double ti = py - i;
sascha@548:         double tj = px - j;
sascha@548: 
sascha@548:         int idx = i*tileWidth + j;
sascha@548: 
sascha@548:         double v1j1 = rasterData[idx];
sascha@548:         double v2j1 = rasterData[idx+1];
sascha@548: 
sascha@548:         idx += tileWidth;
sascha@548: 
sascha@548:         double v1j2 = rasterData[idx];
sascha@548:         double v2j2 = rasterData[idx+1];
sascha@548: 
sascha@548:         double v1 = (1d-tj)*v1j1 + tj*v2j1;
sascha@548:         double v2 = (1d-tj)*v1j2 + tj*v2j2;
sascha@548: 
sascha@548:         return (1d-ti)*v1 + ti*v2;
sascha@548:     }
sascha@548: 
sascha@550:     protected PolynomialSplineFunction [] createSplines() {
sascha@550:         PolynomialSplineFunction [] splines =
sascha@550:             new PolynomialSplineFunction[tileHeight];
sascha@550: 
sascha@550:         double [] x = new double[tileWidth];
sascha@550:         double [] y = new double[tileWidth];
sascha@550: 
sascha@550:         for (int i = 0; i < tileWidth; ++i) {
sascha@550:             x[i] = i + 0.5d;
sascha@550:         }
sascha@550: 
sascha@550:         columns = new double[tileHeight];
sascha@550: 
sascha@550:         SplineInterpolator interpolator = new SplineInterpolator();
sascha@550: 
sascha@550:         for (int i = 0; i < tileHeight; ++i) {
sascha@550:             int row = i*tileWidth;
sascha@550:             for (int j = 0; j < tileWidth; ++j) {
sascha@550:                 y[j] = rasterData[row + j];
sascha@550:             }
sascha@550:             splines[i] = interpolator.interpolate(x, y);
sascha@550:             columns[i] = i + 0.5d;
sascha@550:         }
sascha@550: 
sascha@550:         return splines;
sascha@550:     }
sascha@550: 
sascha@550:     protected synchronized PolynomialSplineFunction [] getSplines() {
sascha@550:         if (splines == null) {
sascha@550:             splines = createSplines();
sascha@550:         }
sascha@550:         return splines;
sascha@550:     }
sascha@550: 
sascha@550:     protected double [] columnValues(double x) {
sascha@550:         PolynomialSplineFunction [] splines = getSplines();
sascha@550:         double [] y = new double[tileHeight];
sascha@550:         try {
sascha@550:             for (int i = 0; i < tileHeight; ++i) {
sascha@550:                 y[i] = splines[i].value(x);
sascha@550:             }
sascha@550:         }
sascha@550:         catch (ArgumentOutsideDomainException aode) {
sascha@550:             log.error(aode.getLocalizedMessage(), aode);
sascha@550:         }
sascha@550:         return y;
sascha@550:     }
sascha@550: 
sascha@550:     private final double evaluateBicubicVertical(double x, double y) {
sascha@550:         try {
sascha@550:             SplineInterpolator interpolator = new SplineInterpolator();
sascha@550:             double [] cs = columnValues(x);
sascha@550:             PolynomialSplineFunction spline = interpolator.interpolate(
sascha@550:                 columns, cs);
sascha@550:             return spline.value(y);
sascha@550:         }
sascha@550:         catch (ArgumentOutsideDomainException aode) {
sascha@550:             log.error(aode.getLocalizedMessage(), aode);
sascha@550:         }
sascha@550:         return 0d;
sascha@550:     }
sascha@550: 
sascha@550:     private final double evaluateBicubicHorizonal(int index, double x) {
sascha@550:         try {
sascha@550:             PolynomialSplineFunction spline = getSplines()[index];
sascha@550:             return spline.value(x);
sascha@550:         }
sascha@550:         catch (ArgumentOutsideDomainException aode) {
sascha@550:             log.error(aode.getLocalizedMessage(), aode);
sascha@550:         }
sascha@550:         return 0d;
sascha@550:     }
sascha@550: 
sascha@550:     public double interpolateBicubic(double px, double py) {
sascha@550: 
sascha@550:         if (px <= 0.5d) { // left
sascha@550:             if (py <= 0.5d) { // upper left
sascha@550:                 return rasterData[0];
sascha@550:             }
sascha@550:             if (py >= tileHeight-0.5d) {
sascha@550:                 return rasterData[rasterData.length-tileWidth];
sascha@550:             }
sascha@550:             return evaluateBicubicVertical(0.5d, py);
sascha@550:         }
sascha@550: 
sascha@550:         if (px >= tileWidth-0.5d) { // right
sascha@550:             if (py < 0.5d) { // upper right edge
sascha@550:                 return rasterData[tileWidth-1];
sascha@550:             }
sascha@550:             if (py > tileHeight-0.5d) { // lower right edge
sascha@550:                 return rasterData[rasterData.length-1];
sascha@550:             }
sascha@550:             return evaluateBicubicVertical(tileHeight-0.5d, py);
sascha@550:         }
sascha@550: 
sascha@550:         if (py <= 0.5d) { // top
sascha@550:             return evaluateBicubicHorizonal(0, px);
sascha@550:         }
sascha@550: 
sascha@550:         if (py >= tileHeight-0.5d) { // bottom
sascha@550:             return evaluateBicubicHorizonal(tileHeight-1, px);
sascha@550:         }
sascha@550: 
sascha@550:         return evaluateBicubicVertical(px, py);
sascha@550:     }
sascha@550: 
sascha@547:     public double getValue(Coordinate coordinate, int interpolationType) {
sascha@547: 
sascha@547:         double px = mx*coordinate.x + bx;
sascha@547:         double py = my*coordinate.y + by;
sascha@547: 
sascha@547:         if (px < 0d || py < 0d || px >= tileWidth || py >= tileHeight) {
sascha@547:             return Double.NaN;
sascha@547:         }
sascha@547: 
sascha@550:         switch (interpolationType) {
sascha@550:             case BILINEAR: return interpolateBilinear(px, py);
sascha@550:             case BICUBIC : return interpolateBicubic (px, py);
sascha@547:         }
sascha@547: 
sascha@547:         int posX = Math.min(tileWidth-1,  (int)Math.round(px));
sascha@547:         int posY = Math.min(tileHeight-1, (int)Math.round(py));
sascha@547: 
sascha@547:         return get(posX, posY);
sascha@547:     }
sascha@549: 
sascha@549:     public static final int getInterpolationType(String key) {
sascha@549:         if (key == null || (key = key.trim().toLowerCase()).length() == 0) {
sascha@549:             return NEAREST_NEIGHBOR;
sascha@549:         }
sascha@549: 
sascha@549:         if ("bilinear".equals(key)) {
sascha@549:             return BILINEAR;
sascha@549:         }
sascha@549: 
sascha@550:         if ("bicubic".equals(key)) {
sascha@550:             return BICUBIC;
sascha@550:         }
sascha@550: 
sascha@549:         return NEAREST_NEIGHBOR;
sascha@549:     }
tim@544: }