bugfix for https://trac.osgeo.org/grass/ticket/475: top bin for FP maps always held just the single max cell (merge from devbr6)

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@37700 15284696-431f-4ddb-bdfa-cd5b030d7da7

raster/r.stats/cell_stats.c

@@ -1,4 +1,5 @@
 #include <stdlib.h>
+#include <grass/glocale.h>
 #include "global.h"
 
 int cell_stats(int fd[], int with_percents, int with_counts,
@@ -47,11 +48,17 @@ int cell_stats(int fd[], int with_percents, int with_counts,
 
 	for (i = 0; i < nfiles; i++) {
 	    if (G_get_c_raster_row(fd[i], cell[i], row) < 0)
-		exit(1);
-	    reset_null_vals(cell[i], ncols);
+		G_fatal_error(_("Unable to read raster map <map %d of %d> row %d"),
+				i+1, nfiles, row);
+
+	    /* include max FP value in nsteps'th bin */
+	    if(is_fp[i])
+		fix_max_fp_val(cell[i], ncols);
+
 	    /* we can't compute hash on null values, so we change all
-	       nulls to max+1, set NULL_CELL to max+1, and later compare with 
-	       NULL_CELL  to chack for nulls */
+	       nulls to max+1, set NULL_CELL to max+1, and later compare
+	       with NULL_CELL to chack for nulls */
+	    reset_null_vals(cell[i], ncols);
 	}
 
 	update_cell_stats(cell, ncols, unit_area);

raster/r.stats/global.h

@@ -24,7 +24,8 @@ int raw_stats(int[], int, int, int);
 int initialize_cell_stats(int);
 int allocate_values(void);
 struct Node *NewNode(double);
-int reset_null_vals(CELL *, int);
+void fix_max_fp_val(CELL *, int);
+void reset_null_vals(CELL *, int);
 int update_cell_stats(CELL **, int, double);
 int node_compare(const void *, const void *);
 int sort_cell_stats(void);

raster/r.stats/main.c

@@ -5,7 +5,7 @@
  * AUTHOR(S):    Michael Shapiro, CERL (original contributor)
  *               Markus Neteler <neteler itc.it>, Brad Douglas <rez touchofmadness.com>,
  *               Huidae Cho <grass4u gmail.com>, Glynn Clements <glynn gclements.plus.com>,
- *               Hamish Bowman <hamish_nospam yahoo.com>,
+ *               Hamish Bowman <hamish_b yahoo.com>,
  *               Jachym Cepicky <jachym les-ejk.cz>, Jan-Oliver Wagner <jan intevation.de>
  * PURPOSE:      calculates the area present in each of the categories of
  *               user-selected raster map layer(s)
@@ -20,9 +20,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include <unistd.h>
-
 #include <grass/glocale.h>
-
 #include "global.h"
 
 char *no_data_str;
@@ -311,10 +309,17 @@ int main(int argc, char *argv[])
 				  name);
 		G_get_fp_range_min_max(&fp_range, &DMIN[nfiles],
 				       &DMAX[nfiles]);
-		G_quant_add_rule(&q, DMIN[nfiles], DMAX[nfiles], 1, nsteps);
+		G_debug(3, "file %2d: dmin=%f  dmax=%f", nfiles, DMIN[nfiles], 
+			DMAX[nfiles]);
+
+		G_quant_add_rule(&q, DMIN[nfiles], DMAX[nfiles], 1, nsteps+1);
+
 		/* set the quant rules for reading the map */
 		G_set_quant_rules(fd[nfiles], &q);
 		G_quant_get_limits(&q, &dmin, &dmax, &min, &max);
+		G_debug(2, "overall: dmin=%f  dmax=%f,  qmin=%d  qmax=%d",
+			dmin, dmax, min, max);
+
 		G_quant_free(&q);
 	    }
 	    else {		/* cats ranges */
@@ -336,6 +341,7 @@ int main(int argc, char *argv[])
 	}
 	else if (NULL_CELL < max + 1)
 	    NULL_CELL = max + 1;
+
 	nfiles++;
     }
 

raster/r.stats/stats.c

@@ -87,16 +87,38 @@ struct Node *NewNode(double area)
     return node;
 }
 
-int reset_null_vals(CELL * cell, int ncols)
+
+/* Essentially, G_quant_add_rule() treats the ranges as half-open,
+ *  i.e. the values range from low (inclusive) to high (exclusive).
+ *  While half-open ranges are a common concept (e.g. floor() behaves
+ *  the same way), the range of a GRASS raster is closed, i.e. both the
+ *  low and high values are inclusive.
+ *  Therefore the quantized max FP cell gets put in the nsteps+1'th bin
+ *  and we need to manually place it back in the previous bin. */
+void fix_max_fp_val(CELL *cell, int ncols)
+{
+    while (ncols-- > 0) {
+	if (cell[ncols] > nsteps)
+	    cell[ncols] = nsteps;
+	    /* { G_debug(5, ". resetting %d to %d\n", cell[ncols], nsteps); } */
+    }
+    return;
+}
+
+
+/* we can't compute hash on null values, so we change all
+ *  nulls to max+1, set NULL_CELL to max+1, and later compare
+ *  with NULL_CELL to chack for nulls */
+void reset_null_vals(CELL *cell, int ncols)
 {
     while (ncols-- > 0) {
 	if (G_is_c_null_value(&cell[ncols]))
 	    cell[ncols] = NULL_CELL;
     }
-
-    return 0;
+    return;
 }
 
+
 int update_cell_stats(CELL ** cell, int ncols, double area)
 {
     register int i;
@@ -108,13 +130,16 @@ int update_cell_stats(CELL ** cell, int ncols, double area)
 	/* copy this cell to an array, compute hash */
 
 	hash = values[0] = cell[0][ncols];
+
 	for (i = 1; i < nfiles; i++)
 	    hash = hash * HASHMOD + (values[i] = cell[i][ncols]);
+
 	if (hash < 0)
 	    hash = -hash;
+
 	hash %= HASHSIZE;
-	/* look it up and update/insert */
 
+	/* look it up and update/insert */
 	if ((q = hashtable[hash]) == NULL) {
 	    hashtable[hash] = NewNode(area);
 	}
@@ -132,6 +157,7 @@ int update_cell_stats(CELL ** cell, int ncols, double area)
 			break;
 		    }
 		}
+
 		if (i == nfiles) {	/* match */
 		    q->count++;
 		    q->area += area;
@@ -192,8 +218,7 @@ int print_node_count(void)
     return 0;
 }
 
-int
-print_cell_stats(char *fmt, int with_percents, int with_counts,
+int print_cell_stats(char *fmt, int with_percents, int with_counts,
 		 int with_areas, int with_labels, char *fs)
 {
     int i, n, nulls_found;