[GRASS-SVN] r67211 - grass/trunk/vector/v.surf.idw

[svn_grass at osgeo.org]

Author: annakrat
Date: 2015-12-17 20:10:41 -0800 (Thu, 17 Dec 2015)
New Revision: 67211

Modified:
   grass/trunk/vector/v.surf.idw/main.c
   grass/trunk/vector/v.surf.idw/proto.h
Log:
v.surf.idw: fix serious bugs - wrong choice of nearest points when using indexing, missing sqrt, see #2820

Modified: grass/trunk/vector/v.surf.idw/main.c
===================================================================
--- grass/trunk/vector/v.surf.idw/main.c	2015-12-18 03:44:10 UTC (rev 67210)
+++ grass/trunk/vector/v.surf.idw/main.c	2015-12-18 04:10:41 UTC (rev 67211)
@@ -62,6 +62,7 @@
     double dist;
     double sum1, sum2, interp_value;
     int n;
+    int max;
     double p;
     struct
     {
@@ -298,6 +299,7 @@
 		else {
 		    pointsfound = 0;
 		    i = 0;
+		    max = 0;
 
 		    if (searchallpoints == 1) {
 			/* If there aren't many sites just check them all to find
@@ -305,7 +307,7 @@
 			for (n = 0; n < ncells; n++)
 			    calculate_distances(shortlistrows[n],
 						shortlistcolumns[n], north,
-						east, &pointsfound);
+						east, &pointsfound, &max);
 		    }
 		    else {
 			radius = 0;
@@ -326,7 +328,7 @@
 					col + search_list[radius]->column;
 				    calculate_distances(searchrow,
 							searchcolumn, north,
-							east, &pointsfound);
+							east, &pointsfound, &max);
 				}
 
 				/* Only if at least one offset is not 0 */
@@ -340,7 +342,7 @@
 					col - search_list[radius]->column;
 				    calculate_distances(searchrow,
 							searchcolumn, north,
-							east, &pointsfound);
+							east, &pointsfound, &max);
 				}
 
 				/* Only if both offsets are not 0 */
@@ -357,7 +359,7 @@
 					calculate_distances(searchrow,
 							    searchcolumn,
 							    north, east,
-							    &pointsfound);
+							    &pointsfound, &max);
 				    }
 				    if (row >= search_list[radius]->row &&
 					col <
@@ -370,7 +372,7 @@
 					calculate_distances(searchrow,
 							    searchcolumn,
 							    north, east,
-							    &pointsfound);
+							    &pointsfound, &max);
 				    }
 				}
 
@@ -386,7 +388,7 @@
 		sum1 = 0.0;
 		sum2 = 0.0;
 		for (n = 0; n < nsearch; n++) {
-		    if ((dist = list[n].dist)) {
+		    if ((dist = sqrt(list[n].dist))) {
 			sum1 += list[n].z / pow(dist, p);
 			sum2 += 1.0 / pow(dist, p);
 		    }
@@ -456,9 +458,9 @@
 }
 
 void calculate_distances(int row, int column, double north,
-			 double east, int *pointsfound)
+			 double east, int *pointsfound, int *max)
 {
-    int j, n, max = 0;
+    int j, n;
     double dx, dy, dist;
     static double maxdist;
 
@@ -474,10 +476,10 @@
 
 	    /* find the maximum distance */
 	    if (i == nsearch) {
-		maxdist = list[max = 0].dist;
+		maxdist = list[*max = 0].dist;
 		for (n = 1; n < nsearch; n++) {
 		    if (maxdist < list[n].dist)
-			maxdist = list[max = n].dist;
+			maxdist = list[*max = n].dist;
 		}
 	    }
 	}
@@ -490,12 +492,12 @@
 
 	    if (dist < maxdist) {
 		/* replace the largest dist */
-		list[max].z = points[row][column][j].z;
-		list[max].dist = dist;
-		maxdist = list[max = 0].dist;
+		list[*max].z = points[row][column][j].z;
+		list[*max].dist = dist;
+		maxdist = list[*max = 0].dist;
 		for (n = 1; n < nsearch; n++) {
 		    if (maxdist < list[n].dist)
-			maxdist = list[max = n].dist;
+			maxdist = list[*max = n].dist;
 		}
 	    }
 

Modified: grass/trunk/vector/v.surf.idw/proto.h
===================================================================
--- grass/trunk/vector/v.surf.idw/proto.h	2015-12-18 03:44:10 UTC (rev 67210)
+++ grass/trunk/vector/v.surf.idw/proto.h	2015-12-18 04:10:41 UTC (rev 67211)
@@ -2,5 +2,5 @@
 void read_sites(const char *, const char *, const char *, int);
 
 void newpoint(double, double, double, int);
-void calculate_distances(int, int, double, double, int *);
+void calculate_distances(int, int, double, double, int *, int *);
 void calculate_distances_noindex(double, double);