# HG changeset patch
# User Sascha L. Teichmann <sascha.teichmann@intevation.de>
# Date 1263256073 0
# Node ID 84ba7cbff791f9ebe7a1e01692379c194b2d90aa
# Parent  0dcf068fb55224233843a7ca715bce2f93d1f8a5
Added bicubic spline interpolation on raster tiles.


geo-backend/trunk@528 c6561f87-3c4e-4783-a992-168aeb5c3f6f

diff -r 0dcf068fb552 -r 84ba7cbff791 geo-backend/ChangeLog
--- a/geo-backend/ChangeLog	Sat Jan 09 16:40:19 2010 +0000
+++ b/geo-backend/ChangeLog	Tue Jan 12 00:27:53 2010 +0000
@@ -1,3 +1,14 @@
+2009-01-12	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
+
+	* src/main/java/de/intevation/gnv/geobackend/sde/datasources/RasterObject.java:
+	  Added bicubic spline interpolation on raster tiles. This
+	  should increase the 'virtual' resolution more than the
+	  bilinear interpolation. To enable it write 'bicubic' into
+	  the configuration file.
+
+	* pom.xml: Added dependency to Apache Common Math 2.0 for
+	  the cubic spline interpolation.
+
 2009-01-09	Sascha L. Teichmann	<sascha.teichmann@intevation.de>
 
 	* src/main/java/de/intevation/gnv/geobackend/sde/datasources/RasterObject.java:
diff -r 0dcf068fb552 -r 84ba7cbff791 geo-backend/pom.xml
--- a/geo-backend/pom.xml	Sat Jan 09 16:40:19 2010 +0000
+++ b/geo-backend/pom.xml	Tue Jan 12 00:27:53 2010 +0000
@@ -62,6 +62,11 @@
       <version>1.5.2</version>
     </dependency>
     <dependency>
+      <groupId>org.apache.commons</groupId>
+      <artifactId>commons-math</artifactId>
+      <version>2.0</version>
+    </dependency>
+    <dependency>
       <groupId>com.vividsolutions</groupId>
       <artifactId>jts</artifactId>
       <version>1.9</version>
diff -r 0dcf068fb552 -r 84ba7cbff791 geo-backend/src/main/java/de/intevation/gnv/geobackend/sde/datasources/RasterObject.java
--- 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;
     }
 }