Add struct FP_stats

Add G_histogram_eq_colors_fp
Change r.colors to use G_histogram_eq_colors_fp() for FP maps


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

include/gis.h

@@ -530,6 +530,15 @@ struct FPRange
     int first_time;		/* whether or not range was updated */
 };
 
+struct FP_stats {
+    int geometric;
+    int flip;
+    int count;
+    DCELL min, max;
+    unsigned long *stats;
+    unsigned long total;
+};
+
 /*
  ** Structure for I/O of 3dview files  (view.c)
  */

include/gisdefs.h

@@ -350,6 +350,8 @@ int G__write_colors(FILE *, struct Colors *);
 /* color_xform.c */
 int G_histogram_eq_colors(struct Colors *, struct Colors *,
 			  struct Cell_stats *);
+void G_histogram_eq_colors_fp(struct Colors *,
+			      struct Colors *, struct FP_stats *);
 int G_log_colors(struct Colors *, struct Colors *, int);
 
 /* commas.c */

lib/gis/color_xform.c

@@ -89,8 +89,6 @@ int G_histogram_eq_colors(struct Colors *dst,
     prev = 0;
     first = 1;
 
-    G_get_d_raster_color(&min, &red, &grn, &blu, src);
-
     G_rewind_cell_stats(statf);
     while (G_next_cell_stat(&cat, &count, statf)) {
 	int red2, grn2, blu2;
@@ -118,6 +116,75 @@ int G_histogram_eq_colors(struct Colors *dst,
 }
 
 /*!
+ * \brief make histogram-stretched version of existing color table (FP version)
+ *
+ * Generates a histogram
+ * contrast-stretched color table that goes from the histogram
+ * information in the FP_stats structure <b>statf.</b>  (See
+ * Raster_Histograms).
+ *
+ *  \param dst
+ *  \param src
+ *  \param statf
+ *  \return void
+ */
+
+void G_histogram_eq_colors_fp(struct Colors *dst,
+			      struct Colors *src, struct FP_stats *statf)
+{
+    DCELL min, max;
+    int red, grn, blu;
+    unsigned long sum;
+    DCELL val;
+    int first;
+    int i;
+
+    G_init_colors(dst);
+
+    G_get_d_color_range(&min, &max, src);
+
+    G_get_default_color(&red, &grn, &blu, src);
+    G_set_default_color(red, grn, blu, dst);
+
+    G_get_null_value_color(&red, &grn, &blu, src);
+    G_set_null_value_color(red, grn, blu, dst);
+
+    if (!statf->total)
+	return;
+
+    sum = 0;
+    first = 1;
+
+    for (i = 0; i <= statf->count; i++) {
+	int red2, grn2, blu2;
+	DCELL val2, x;
+
+	val2 = statf->min + (statf->max - statf->min) * i / statf->count;
+	if (statf->geometric)
+	    val2 = exp(val2);
+	if (statf->flip)
+	    val2 = -val2;
+
+	x = min + (max - min) * sum / statf->total;
+	G_get_d_raster_color(&x, &red2, &grn2, &blu2, src);
+
+	if (!first)
+	    G_add_d_raster_color_rule(&val, red, grn, blu, &val2, red2, grn2, blu2, dst);
+	first = 0;
+
+	if (i == statf->count)
+	    break;
+
+	sum += statf->stats[i];
+
+	val = val2;
+	red = red2;
+	grn = grn2;
+	blu = blu2;
+    }
+}
+
+/*!
  * \brief make logarithmically-scaled version of an existing color table
  *
  *  \param dst

raster/r.colors/local_proto.h

@@ -23,7 +23,10 @@
 #include <grass/gis.h>
 
 /* stats.c */
-int get_stats(char *, char *, struct Cell_stats *);
+int get_stats(const char *, const char *, struct Cell_stats *);
+void get_fp_stats(const char *name, const char *mapset,
+		  struct FP_stats *statf,
+		  DCELL min, DCELL max, int geometric);
 
 /* main.c */
 int main(int, char *[]);

raster/r.colors/main.c

@@ -342,9 +342,16 @@ int main(int argc, char **argv)
 	G_invert_colors(&colors);
 
     if (flag.e->answer) {
-	if (!have_stats)
-	    have_stats = get_stats(name, mapset, &statf);
-	G_histogram_eq_colors(&colors_tmp, &colors, &statf);
+	if (fp) {
+	    struct FP_stats fpstats;
+	    get_fp_stats(name, mapset, &fpstats, min, max, flag.g->answer);
+	    G_histogram_eq_colors_fp(&colors_tmp, &colors, &fpstats);
+	}
+	else {
+	    if (!have_stats) 
+		have_stats = get_stats(name, mapset, &statf);
+	    G_histogram_eq_colors(&colors_tmp, &colors, &statf);
+	}
 	colors = colors_tmp;
     }
 

raster/r.colors/stats.c

@@ -17,11 +17,13 @@
  *
  ***************************************************************************/
 
+#include <math.h>
 #include <stdlib.h>
 #include <grass/gis.h>
 #include <grass/glocale.h>
+#include "local_proto.h"
 
-int get_stats(char *name, char *mapset, struct Cell_stats *statf)
+int get_stats(const char *name, const char *mapset, struct Cell_stats *statf)
 {
     CELL *cell;
     int row, nrows, ncols;
@@ -48,3 +50,72 @@ int get_stats(char *name, char *mapset, struct Cell_stats *statf)
 
     return 1;
 }
+
+void get_fp_stats(const char *name, const char *mapset,
+		  struct FP_stats *statf,
+		  DCELL min, DCELL max, int geometric)
+{
+    DCELL *dcell;
+    int row, col, nrows, ncols;
+    int fd;
+
+    if ((fd = G_open_cell_old(name, mapset)) < 0)
+	G_fatal_error("error opening map <%s@%s>", name, mapset);
+
+    dcell = G_allocate_d_raster_buf();
+    nrows = G_window_rows();
+    ncols = G_window_cols();
+
+    statf->geometric = geometric;
+    statf->flip = 0;
+
+    if (statf->geometric) {
+	if (min * max < 0)
+	    G_fatal_error(_("Cannot use logarithmic scaling if range includes zero"));
+
+	if (min < 0) {
+	    statf->flip = 1;
+	    min = -min;
+	    max = -max;
+	}
+
+	min = log(min);
+	max = log(max);
+    }
+
+    statf->count = 1000;
+    statf->min = min;
+    statf->max = max;
+    statf->stats = G_calloc(statf->count, sizeof(unsigned long));
+    statf->total = 0;
+
+    G_message(_("Reading %s ..."), name);
+    for (row = 0; row < nrows; row++) {
+	G_percent(row, nrows, 2);
+
+	if (G_get_d_raster_row(fd, dcell, row) < 0)
+	    G_fatal_error("error reading map <%s@%s>", name, mapset);
+
+	for (col = 0; col < ncols; col++) {
+	    DCELL x;
+	    int i;
+
+	    if (G_is_d_null_value(&dcell[col]))
+		continue;
+
+	    x = dcell[col];
+	    if (statf->flip)
+		x = -x;
+	    if (statf->geometric)
+		x = log(x);
+
+	    i = (int) floor(statf->count * (x - statf->min) / (statf->max - statf->min));
+	    statf->stats[i]++;
+	    statf->total++;
+	}
+    }
+
+    G_percent(row, nrows, 2);
+    G_close_cell(fd);
+    G_free(dcell);
+}