Mercurial > dive4elements > gnv-client
view geo-backend/src/main/java/de/intevation/gnv/geobackend/sde/datasources/RasterObject.java @ 1127:ebeb56428409
Added license headers and license file.
geo-backend/trunk@1261 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Ingo Weinzierl <ingo.weinzierl@intevation.de> |
---|---|
date | Tue, 02 Nov 2010 17:52:22 +0000 |
parents | b757def3ff55 |
children |
line wrap: on
line source
/* * Copyright (c) 2010 by Intevation GmbH * * This program is free software under the LGPL (>=v2.1) * Read the file LGPL.txt coming with the software for details * or visit http://www.gnu.org/licenses/ if it does not exist. */ package de.intevation.gnv.geobackend.sde.datasources; import com.vividsolutions.jts.geom.Coordinate; import java.io.Serializable; import org.apache.commons.math.ArgumentOutsideDomainException; import org.apache.commons.math.analysis.interpolation.SplineInterpolator; import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; import org.apache.log4j.Logger; /** * @author <a href="mailto:tim.englich@intevation.de">Tim Englich</a> * @author <a href="mailto:sascha.teichmann@intevation.de">Sascha L. Teichmann</a> */ public class RasterObject implements Serializable { private static Logger log = Logger.getLogger(RasterObject.class); 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; private double my; private double by; private int columnIndex; private int rowIndex; private int tileWidth; private int tileHeight; private double [] rasterData; private transient PolynomialSplineFunction [] splines; private transient double [] columns; public RasterObject() { } public RasterObject( double mx, double bx, double my, double by, int columnIndex, int rowIndex, double [] rasterData, int tileWidth, int tileHeight ){ this.mx = mx; this.bx = bx; this.my = my; this.by = by; this.columnIndex = columnIndex; this.rowIndex = rowIndex; this.rasterData = rasterData; this.tileWidth = tileWidth; this.tileHeight = tileHeight; } public boolean contains(Coordinate coordinate) { double px = mx*coordinate.x + bx; double py = my*coordinate.y + by; return px >= 0d && py >= 0d && px < tileWidth && py < tileHeight; } public final double get(int posX, int posY) { return rasterData[posY*tileWidth + posX]; } public int getTileWidth() { return tileWidth; } public int getTileHeight() { return tileWidth; } public int getColumnIndex() { return columnIndex; } public int getRowIndex() { return rowIndex; } public double getValue(Coordinate coordinate) { return getValue(coordinate, NEAREST_NEIGHBOR); } public final double interpolateBilinear(double px, double py) { if (px < 0.5d) { // left border if (py < 0.5d) { // upper left edge return rasterData[0]; } if (py > tileHeight-0.5d) { // lower left edge return rasterData[rasterData.length-tileWidth]; } // center left int i = (int)(py -= 0.5d); double t = py - i; i *= tileWidth; return (1d-t)*rasterData[i] + t*rasterData[i+tileWidth]; } if (px > tileWidth-0.5d) { // right border if (py < 0.5d) { // upper right edge return rasterData[tileWidth-1]; } if (py > tileHeight-0.5d) { // lower right edge return rasterData[rasterData.length-1]; } // center right int i = (int)(py -= 0.5d); double t = py - i; i = i*tileWidth + tileWidth-1; return (1d-t)*rasterData[i] + t*rasterData[i+tileWidth]; } if (py < 0.5d) { // top border int i = (int)(px -= 0.5d); double t = px - i; return (1d-t)*rasterData[i] + t*rasterData[i+1]; } if (py > tileHeight-0.5d) { // bottom border int i = (int)(px -= 0.5d); double t = px - i; i += rasterData.length-tileWidth; return (1d-t)*rasterData[i] + t*rasterData[i+1]; } // center int i = (int)(py -= 0.5d); int j = (int)(px -= 0.5d); double ti = py - i; double tj = px - j; int idx = i*tileWidth + j; double v1j1 = rasterData[idx]; double v2j1 = rasterData[idx+1]; idx += tileWidth; double v1j2 = rasterData[idx]; double v2j2 = rasterData[idx+1]; double v1 = (1d-tj)*v1j1 + tj*v2j1; double v2 = (1d-tj)*v1j2 + tj*v2j2; 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; double py = my*coordinate.y + by; if (px < 0d || py < 0d || px >= tileWidth || py >= tileHeight) { return Double.NaN; } switch (interpolationType) { case BILINEAR: return interpolateBilinear(px, py); case BICUBIC : return interpolateBicubic (px, py); } int posX = Math.min(tileWidth-1, (int)Math.round(px)); int posY = Math.min(tileHeight-1, (int)Math.round(py)); return get(posX, posY); } public static final int getInterpolationType(String key) { if (key == null || (key = key.trim().toLowerCase()).length() == 0) { return NEAREST_NEIGHBOR; } if ("bilinear".equals(key)) { return BILINEAR; } if ("bicubic".equals(key)) { return BICUBIC; } return NEAREST_NEIGHBOR; } }