[GRASS-SVN] r44058 - grass/trunk/lib/raster

[svn_grass at osgeo.org]

Author: mmetz
Date: 2010-10-27 23:36:13 -0700 (Wed, 27 Oct 2010)
New Revision: 44058

Modified:
   grass/trunk/lib/raster/interp.c
Log:
add lanczos and cubic bspline interpolation to rasterlib

Modified: grass/trunk/lib/raster/interp.c
===================================================================
--- grass/trunk/lib/raster/interp.c	2010-10-28 02:53:30 UTC (rev 44057)
+++ grass/trunk/lib/raster/interp.c	2010-10-28 06:36:13 UTC (rev 44058)
@@ -11,9 +11,12 @@
  * \author Original author CERL
  */
 
+#include <math.h>
 #include <grass/gis.h>
 #include <grass/raster.h>
 
+#define LANCZOS_FILTER(x) ((2 * sin(x) * sin((x) / 2)) / ((x) * (x)))
+
 DCELL Rast_interp_linear(double u, DCELL c0, DCELL c1)
 {
     return u * (c1 - c0) + c0;
@@ -48,3 +51,58 @@
 
     return Rast_interp_cubic(v, c0, c1, c2, c3);
 }
+
+DCELL Rast_interp_lanczos(double u, double v, DCELL *c)
+{
+    int i, j;
+    double uweight, vweight[5], d;
+    DCELL result = 0;
+
+    for (i = 0; i < 5; i++) {
+	d = u - i + 2;
+	if (d == 0)
+	    uweight = 1;
+	else {
+	    d *= M_PI;
+	    uweight = LANCZOS_FILTER(d);
+	}
+
+	for (j = 0; j < 5; j++) {
+	    if (i == 0) {
+		d = v - j + 2;
+		if (d == 0)
+		    vweight[j] = 1;
+		else {
+		    d *= M_PI;
+		    vweight[j] = LANCZOS_FILTER(d);
+		}
+	    }
+
+	    result += *(c++) * uweight * vweight[j];
+	}
+    }
+    
+    return result;
+}
+
+DCELL Rast_interp_cubic_bspline(double u, DCELL c0, DCELL c1, DCELL c2, DCELL c3)
+{
+    return (u * (u * (u * ( 1 * c3 - 3 * c2 + 3 * c1 - 1 * c0) +
+			  ( 0 * c3 + 3 * c2 - 6 * c1 + 3 * c0)) +
+			  ( 0 * c3 + 3 * c2 + 0 * c1 - 3 * c0)) +
+			    0 * c3 + 1 * c2 + 4 * c1 + 1 * c0) / 6;
+}
+
+DCELL Rast_interp_bicubic_bspline(double u, double v,
+			  DCELL c00, DCELL c01, DCELL c02, DCELL c03,
+			  DCELL c10, DCELL c11, DCELL c12, DCELL c13,
+			  DCELL c20, DCELL c21, DCELL c22, DCELL c23,
+			  DCELL c30, DCELL c31, DCELL c32, DCELL c33)
+{
+    DCELL c0 = Rast_interp_bspline_cubic(u, c00, c01, c02, c03);
+    DCELL c1 = Rast_interp_bspline_cubic(u, c10, c11, c12, c13);
+    DCELL c2 = Rast_interp_bspline_cubic(u, c20, c21, c22, c23);
+    DCELL c3 = Rast_interp_bspline_cubic(u, c30, c31, c32, c33);
+
+    return Rast_interp_bspline_cubic(v, c0, c1, c2, c3);
+}