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

[svn_grass at osgeo.org]

Author: mmetz
Date: 2010-10-29 09:02:20 -0700 (Fri, 29 Oct 2010)
New Revision: 44082

Modified:
   grass/trunk/lib/raster/interp.c
Log:
avoid fp calculation errors as far as possible

Modified: grass/trunk/lib/raster/interp.c
===================================================================
--- grass/trunk/lib/raster/interp.c	2010-10-29 13:09:35 UTC (rev 44081)
+++ grass/trunk/lib/raster/interp.c	2010-10-29 16:02:20 UTC (rev 44082)
@@ -55,51 +55,49 @@
     double uweight[5], vweight[5], d, d_pi;
     double sind, sincd1, sincd2;
 
-    d = v + 2;
-    d_pi = d * M_PI;
+    d_pi = v * M_PI;
     sind = 2 * sin(d_pi);
     sincd1 = sind * sin(d_pi / 2);
-    vweight[0] = (d == 0 ? 1 : sincd1 / (d_pi * d_pi));
+    vweight[2] = (v == 0 ? 1 : sincd1 / (d_pi * d_pi));
 
-    d -= 1.;
+    d = v + 2;
     d_pi = d * M_PI;
-    sincd2 = -sind * sin(d_pi / 2);
-    vweight[1] = (d == 0 ? 1 : sincd2 / (d_pi * d_pi));
+    vweight[0] = (d == 0 ? 1 : -sincd1 / (d_pi * d_pi));
 
-    d -= 1.;
+    d = v + 1.;
     d_pi = d * M_PI;
-    vweight[2] = (d == 0 ? 1 : -sincd1 / (d_pi * d_pi));
+    sincd2 = sind * sin(d_pi / 2);
+    vweight[1] = (d == 0 ? 1 : -sincd2 / (d_pi * d_pi));
 
-    d -= 1.;
+    d = v - 1.;
     d_pi = d * M_PI;
-    vweight[3] = (d == 0 ? 1 : -sincd2 / (d_pi * d_pi));
+    vweight[3] = (d == 0 ? 1 : sincd2 / (d_pi * d_pi));
 
-    d -= 1.;
+    d = v - 2.;
     d_pi = d * M_PI;
-    vweight[4] = (d == 0 ? 1 : sincd1 / (d_pi * d_pi));
+    vweight[4] = (d == 0 ? 1 : -sincd1 / (d_pi * d_pi));
 
-    d = u + 2;
-    d_pi = d * M_PI;
+    d_pi = u * M_PI;
     sind = 2 * sin(d_pi);
     sincd1 = sind * sin(d_pi / 2);
-    uweight[0] = (d == 0 ? 1 : sincd1 / (d_pi * d_pi));
+    uweight[2] = (u == 0 ? 1 : sincd1 / (d_pi * d_pi));
 
-    d -= 1.;
+    d = u + 2;
     d_pi = d * M_PI;
-    sincd2 = -sind * sin(d_pi / 2);
-    uweight[1] = (d == 0 ? 1 : sincd2 / (d_pi * d_pi));
+    uweight[0] = (d == 0 ? 1 : -sincd1 / (d_pi * d_pi));
 
-    d -= 1.;
+    d = u + 1.;
     d_pi = d * M_PI;
-    uweight[2] = (d == 0 ? 1 : -sincd1 / (d_pi * d_pi));
+    sincd2 = sind * sin(d_pi / 2);
+    uweight[1] = (d == 0 ? 1 : -sincd2 / (d_pi * d_pi));
 
-    d -= 1.;
+    d = u - 1.;
     d_pi = d * M_PI;
-    uweight[3] = (d == 0 ? 1 : -sincd2 / (d_pi * d_pi));
+    uweight[3] = (d == 0 ? 1 : sincd2 / (d_pi * d_pi));
 
-    d -= 1.;
+    d = u - 2.;
     d_pi = d * M_PI;
-    uweight[4] = (d == 0 ? 1 : sincd1 / (d_pi * d_pi));
+    uweight[4] = (d == 0 ? 1 : -sincd1 / (d_pi * d_pi));
 
     return ((c[0] * vweight[0] + c[1] * vweight[1] + c[2] * vweight[2]
 	    + c[3] * vweight[3] + c[4] * vweight[4]) * uweight[0] +