[GRASS-SVN] r70879 - grass-addons/grass7/raster/r.fill.gaps
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Apr 14 15:31:39 PDT 2017
Author: wenzeslaus
Date: 2017-04-14 15:31:39 -0700 (Fri, 14 Apr 2017)
New Revision: 70879
Modified:
grass-addons/grass7/raster/r.fill.gaps/cell_funcs.c
grass-addons/grass7/raster/r.fill.gaps/cell_funcs.h
grass-addons/grass7/raster/r.fill.gaps/main.c
Log:
r.fill.gaps: run indent script, remove tabs from first comment
Modified: grass-addons/grass7/raster/r.fill.gaps/cell_funcs.c
===================================================================
--- grass-addons/grass7/raster/r.fill.gaps/cell_funcs.c 2017-04-14 22:25:31 UTC (rev 70878)
+++ grass-addons/grass7/raster/r.fill.gaps/cell_funcs.c 2017-04-14 22:31:39 UTC (rev 70879)
@@ -1,18 +1,20 @@
+
/***************************************************************************
*
- * MODULE: r.fill.gaps (commandline)
- * FILE: cell_funcs.c
- * AUTHOR(S): Benjamin Ducke
- * PURPOSE: Provide cell type dependent functions for cell data operations.
- * At program initialization, a function pointer is set to point to
- * the right one of the three alternatives.
- * This function pointer can later be used to work with cell data of
- * different types without having to go through any switching logics.
+ * MODULE: r.fill.gaps (commandline)
+ * FILE: cell_funcs.c
+ * AUTHOR(S): Benjamin Ducke
+ * PURPOSE: Provide cell type dependent functions for cell data operations.
+ * At program initialization, a function pointer is set to point to
+ * the right one of the three alternatives.
+ * This function pointer can later be used to work with cell data of
+ * different types without having to go through any switching logics.
*
- * COPYRIGHT: (C) 2011 by the GRASS Development Team
+ * COPYRIGHT: (C) 2011 by the GRASS Development Team
*
- * This program is free software under the GPL (>=v2)
- * Read the file COPYING that comes with GRASS for details.
+ * This program is free software under the GPL (>=v2)
+ * Read the file COPYING that comes with GRASS for details.
+ *
****************************************************************************
*/
@@ -28,16 +30,22 @@
/*
* Write cell values.
*/
-void write_cell_value_c ( void *cell_output, void *cell_input ) {
- G_set_raster_value_c(cell_output,G_get_raster_value_c(cell_input,IN_TYPE), OUT_TYPE);
+void write_cell_value_c(void *cell_output, void *cell_input)
+{
+ G_set_raster_value_c(cell_output,
+ G_get_raster_value_c(cell_input, IN_TYPE), OUT_TYPE);
}
-void write_cell_value_f ( void *cell_output, void *cell_input ) {
- G_set_raster_value_f(cell_output,G_get_raster_value_f(cell_input,IN_TYPE), OUT_TYPE);
+void write_cell_value_f(void *cell_output, void *cell_input)
+{
+ G_set_raster_value_f(cell_output,
+ G_get_raster_value_f(cell_input, IN_TYPE), OUT_TYPE);
}
-void write_cell_value_d ( void *cell_output, void *cell_input ) {
- G_set_raster_value_d(cell_output,G_get_raster_value_d(cell_input,IN_TYPE), OUT_TYPE);
+void write_cell_value_d(void *cell_output, void *cell_input)
+{
+ G_set_raster_value_d(cell_output,
+ G_get_raster_value_d(cell_input, IN_TYPE), OUT_TYPE);
}
@@ -45,16 +53,19 @@
* Write a double value into a cell (truncates for CELL type output).
*/
-void write_double_value_c ( void *cell, double val ) {
- G_set_raster_value_c(cell,(CELL)val, OUT_TYPE);
+void write_double_value_c(void *cell, double val)
+{
+ G_set_raster_value_c(cell, (CELL) val, OUT_TYPE);
}
-void write_double_value_f ( void *cell, double val ) {
- G_set_raster_value_f(cell,(FCELL)val, OUT_TYPE);
+void write_double_value_f(void *cell, double val)
+{
+ G_set_raster_value_f(cell, (FCELL) val, OUT_TYPE);
}
-void write_double_value_d ( void *cell, double val ) {
- G_set_raster_value_d(cell,(DCELL)val, OUT_TYPE);
+void write_double_value_d(void *cell, double val)
+{
+ G_set_raster_value_d(cell, (DCELL) val, OUT_TYPE);
}
@@ -62,32 +73,38 @@
/*
* Check for "no data".
*/
-int is_null_value_c ( void *cell ) {
- return (G_is_c_null_value ((CELL*)cell));
+int is_null_value_c(void *cell)
+{
+ return (G_is_c_null_value((CELL *) cell));
}
-int is_null_value_f ( void *cell ) {
- return (G_is_f_null_value ((FCELL*)cell));
+int is_null_value_f(void *cell)
+{
+ return (G_is_f_null_value((FCELL *) cell));
}
-int is_null_value_d ( void *cell ) {
- return (G_is_d_null_value ((DCELL*)cell));
+int is_null_value_d(void *cell)
+{
+ return (G_is_d_null_value((DCELL *) cell));
}
/*
* Set consecutive cells to "no data".
*/
-void set_null_c (void* cells, unsigned long count) {
- G_set_c_null_value ((CELL*)cells, count);
+void set_null_c(void *cells, unsigned long count)
+{
+ G_set_c_null_value((CELL *) cells, count);
}
-void set_null_f (void* cells, unsigned long count) {
- G_set_f_null_value ((FCELL*)cells, count);
+void set_null_f(void *cells, unsigned long count)
+{
+ G_set_f_null_value((FCELL *) cells, count);
}
-void set_null_d (void* cells, unsigned long count) {
- G_set_d_null_value ((DCELL*)cells, count);
+void set_null_d(void *cells, unsigned long count)
+{
+ G_set_d_null_value((DCELL *) cells, count);
}
@@ -97,44 +114,45 @@
* I.e. right after IN_TYPE and OUT_TYPE have been set
* to valid values;
*/
-void init_cell_funcs () {
+void init_cell_funcs()
+{
- CELL_ERR_SIZE = sizeof (FCELL);
+ CELL_ERR_SIZE = sizeof(FCELL);
- if ( IN_TYPE == CELL_TYPE ) {
- CELL_IN_SIZE = sizeof (CELL);
- CELL_IN_PTR_SIZE = sizeof (CELL*);
- IS_NULL = &is_null_value_c;
- }
- if ( IN_TYPE == FCELL_TYPE ) {
- CELL_IN_SIZE = sizeof (FCELL);
- CELL_IN_PTR_SIZE = sizeof (FCELL*);
- IS_NULL = &is_null_value_f;
- }
- if ( IN_TYPE == DCELL_TYPE ) {
- CELL_IN_SIZE = sizeof (DCELL);
- CELL_IN_PTR_SIZE = sizeof (DCELL*);
- IS_NULL = &is_null_value_d;
- }
- if ( OUT_TYPE == CELL_TYPE ) {
- CELL_OUT_SIZE = sizeof (CELL);
- CELL_OUT_PTR_SIZE = sizeof (CELL*);
- WRITE_CELL_VAL = &write_cell_value_c;
- WRITE_DOUBLE_VAL = &write_double_value_c;
- SET_NULL = &set_null_c;
- }
- if ( OUT_TYPE == FCELL_TYPE ) {
- CELL_OUT_SIZE = sizeof (FCELL);
- CELL_OUT_PTR_SIZE = sizeof (FCELL*);
- WRITE_CELL_VAL = &write_cell_value_f;
- WRITE_DOUBLE_VAL = &write_double_value_f;
- SET_NULL = &set_null_f;
- }
- if ( OUT_TYPE == DCELL_TYPE ) {
- CELL_OUT_SIZE = sizeof (DCELL);
- CELL_OUT_PTR_SIZE = sizeof (DCELL*);
- WRITE_CELL_VAL = &write_cell_value_d;
- WRITE_DOUBLE_VAL = &write_double_value_d;
- SET_NULL = &set_null_d;
- }
+ if (IN_TYPE == CELL_TYPE) {
+ CELL_IN_SIZE = sizeof(CELL);
+ CELL_IN_PTR_SIZE = sizeof(CELL *);
+ IS_NULL = &is_null_value_c;
+ }
+ if (IN_TYPE == FCELL_TYPE) {
+ CELL_IN_SIZE = sizeof(FCELL);
+ CELL_IN_PTR_SIZE = sizeof(FCELL *);
+ IS_NULL = &is_null_value_f;
+ }
+ if (IN_TYPE == DCELL_TYPE) {
+ CELL_IN_SIZE = sizeof(DCELL);
+ CELL_IN_PTR_SIZE = sizeof(DCELL *);
+ IS_NULL = &is_null_value_d;
+ }
+ if (OUT_TYPE == CELL_TYPE) {
+ CELL_OUT_SIZE = sizeof(CELL);
+ CELL_OUT_PTR_SIZE = sizeof(CELL *);
+ WRITE_CELL_VAL = &write_cell_value_c;
+ WRITE_DOUBLE_VAL = &write_double_value_c;
+ SET_NULL = &set_null_c;
+ }
+ if (OUT_TYPE == FCELL_TYPE) {
+ CELL_OUT_SIZE = sizeof(FCELL);
+ CELL_OUT_PTR_SIZE = sizeof(FCELL *);
+ WRITE_CELL_VAL = &write_cell_value_f;
+ WRITE_DOUBLE_VAL = &write_double_value_f;
+ SET_NULL = &set_null_f;
+ }
+ if (OUT_TYPE == DCELL_TYPE) {
+ CELL_OUT_SIZE = sizeof(DCELL);
+ CELL_OUT_PTR_SIZE = sizeof(DCELL *);
+ WRITE_CELL_VAL = &write_cell_value_d;
+ WRITE_DOUBLE_VAL = &write_double_value_d;
+ SET_NULL = &set_null_d;
+ }
}
Modified: grass-addons/grass7/raster/r.fill.gaps/cell_funcs.h
===================================================================
--- grass-addons/grass7/raster/r.fill.gaps/cell_funcs.h 2017-04-14 22:25:31 UTC (rev 70878)
+++ grass-addons/grass7/raster/r.fill.gaps/cell_funcs.h 2017-04-14 22:31:39 UTC (rev 70879)
@@ -1,22 +1,25 @@
+
/***************************************************************************
*
- * MODULE: r.fill.gaps (commandline)
- * FILE: cell_funcs.h
- * AUTHOR(S): Benjamin Ducke
- * PURPOSE: See cell_funcs.c
+ * MODULE: r.fill.gaps (commandline)
+ * FILE: cell_funcs.h
+ * AUTHOR(S): Benjamin Ducke
+ * PURPOSE: See cell_funcs.c
*
- * COPYRIGHT: (C) 2011 by the GRASS Development Team
+ * COPYRIGHT: (C) 2011 by the GRASS Development Team
*
- * This program is free software under the GPL (>=v2)
- * Read the file COPYING that comes with GRASS for details.
+ * This program is free software under the GPL (>=v2)
+ * Read the file COPYING that comes with GRASS for details.
+ *
****************************************************************************
*/
#ifndef CELL_FUNCS_H
#define CELL_FUNCS_H
-RASTER_MAP_TYPE IN_TYPE; /* stuff for cell input and output data */
+RASTER_MAP_TYPE IN_TYPE; /* stuff for cell input and output data */
RASTER_MAP_TYPE OUT_TYPE;
+
/* below are the sizes (in bytes) for the different GRASS cell types */
unsigned char CELL_IN_SIZE;
unsigned char CELL_IN_PTR_SIZE;
@@ -24,9 +27,9 @@
unsigned char CELL_OUT_PTR_SIZE;
unsigned char CELL_ERR_SIZE;
-void (*WRITE_CELL_VAL) (void*, void*); /* write a cell value of any type into an output cell */
-void (*WRITE_DOUBLE_VAL) (void*, double); /* write a double value into an output cell */
-int (*IS_NULL) (void*); /* check if a cell is "no data" */
-void (*SET_NULL) (void*, unsigned long); /* write null value(s) */
+void (*WRITE_CELL_VAL) (void *, void *); /* write a cell value of any type into an output cell */
+void (*WRITE_DOUBLE_VAL) (void *, double); /* write a double value into an output cell */
+int (*IS_NULL) (void *); /* check if a cell is "no data" */
+void (*SET_NULL) (void *, unsigned long); /* write null value(s) */
#endif /* CELL_FUNCS_H */
Modified: grass-addons/grass7/raster/r.fill.gaps/main.c
===================================================================
--- grass-addons/grass7/raster/r.fill.gaps/main.c 2017-04-14 22:25:31 UTC (rev 70878)
+++ grass-addons/grass7/raster/r.fill.gaps/main.c 2017-04-14 22:31:39 UTC (rev 70879)
@@ -1,43 +1,44 @@
/***************************************************************************
*
- * MODULE: r.fill.gaps (commandline)
- * AUTHOR(S): Benjamin Ducke
- * PURPOSE: Rapidly fill "no data" cells in a raster map by simple interpolation.
+ * MODULE: r.fill.gaps (commandline)
+ * AUTHOR(S): Benjamin Ducke
+ * PURPOSE: Rapidly fill "no data" cells in a raster map by simple interpolation.
*
- * COPYRIGHT: (C) 2011-2015 by the GRASS Development Team
+ * COPYRIGHT: (C) 2011-2015 by the GRASS Development Team
*
- * This program is free software under the GPL (>=v2)
- * Read the file COPYING that comes with GRASS for details.
+ * This program is free software under the GPL (>=v2)
+ * Read the file COPYING that comes with GRASS for details.
+ *
****************************************************************************
*/
/*
-DOCS:
-- docs need a lot of updates!
-- dropped medoid (since we always have an odd number of cells, median=medoid).
-- flag -p(reserve) reversed (preservation is now default)
-- flag for including center cell has been dropped; center always included
-- "method" renamed to "mode"
-- by default, neighborhood size is now assumed to be in cells,
- unless -m(ap units) is given
-- only -m(ap units) will result in an exact, circular neighborhood;
- if (default) cells size specification is used, then the neighborhood
- shape is a less exact rectangle; in that case the neighborhood dimensions
- will also be different for x and y if the cell dimensions are different
-- lat/lon data is allowed, but distance measure for -m(ap units) is straight line!
-- cell weights are now normalized to be in range [0;1]
-- center cell weight for "wmean" is now "1.0"
+ DOCS:
+ - docs need a lot of updates!
+ - dropped medoid (since we always have an odd number of cells, median=medoid).
+ - flag -p(reserve) reversed (preservation is now default)
+ - flag for including center cell has been dropped; center always included
+ - "method" renamed to "mode"
+ - by default, neighborhood size is now assumed to be in cells,
+ unless -m(ap units) is given
+ - only -m(ap units) will result in an exact, circular neighborhood;
+ if (default) cells size specification is used, then the neighborhood
+ shape is a less exact rectangle; in that case the neighborhood dimensions
+ will also be different for x and y if the cell dimensions are different
+ - lat/lon data is allowed, but distance measure for -m(ap units) is straight line!
+ - cell weights are now normalized to be in range [0;1]
+ - center cell weight for "wmean" is now "1.0"
-BUGS:
-- mode, median and mode produce large areas of "no data" where there is input data!!!
+ BUGS:
+ - mode, median and mode produce large areas of "no data" where there is input data!!!
-NEXT VERSION:
-- add lat/lon distance and cost map distance measures (lat/lon currently throws an error).
-- allow user to set which cell value should be filled (i.e. interpreted as "no data")
+ NEXT VERSION:
+ - add lat/lon distance and cost map distance measures (lat/lon currently throws an error).
+ - allow user to set which cell value should be filled (i.e. interpreted as "no data")
-*/
+ */
#include <math.h>
#include <string.h>
@@ -51,8 +52,8 @@
/* Global variables :-) */
-double **WEIGHTS; /* neighborhood weights */
-double SUM_WEIGHTS; /* sum of weights in all cells (with or without data) of neighborhood */
+double **WEIGHTS; /* neighborhood weights */
+double SUM_WEIGHTS; /* sum of weights in all cells (with or without data) of neighborhood */
unsigned long WINDOW_WIDTH = 0; /* dimensions of neighborhood window */
unsigned long WINDOW_HEIGHT = 0;
@@ -69,19 +70,20 @@
/* holds statistics of cells within the neighborhood */
typedef struct
{
- unsigned long num_values; /* number of cells with values in input raster */
- double *values; /* individual values of all cells */
- double *weights; /* individual weights of all cells */
- double result; /* weighted mean of values in entire neighborhood */
- double certainty; /* certainty measure, always between 0 (lowest) and 1 (highest) */
- unsigned long *frequencies; /* frequency count for each value */
- double *overwrite; /* will be set to non-null to overwrite the statistical result with this value */
+ unsigned long num_values; /* number of cells with values in input raster */
+ double *values; /* individual values of all cells */
+ double *weights; /* individual weights of all cells */
+ double result; /* weighted mean of values in entire neighborhood */
+ double certainty; /* certainty measure, always between 0 (lowest) and 1 (highest) */
+ unsigned long *frequencies; /* frequency count for each value */
+ double *overwrite; /* will be set to non-null to overwrite the statistical result with this value */
} stats_struct;
/* function pointers for operation modes */
-void (*GET_STATS) (unsigned long, unsigned long, double, double, int, stats_struct*);
-void (*COLLECT_DATA) (double, double, double, double, stats_struct*);
+void (*GET_STATS) (unsigned long, unsigned long, double, double, int,
+ stats_struct *);
+void (*COLLECT_DATA) (double, double, double, double, stats_struct *);
/*
@@ -89,49 +91,49 @@
* (in bytes) that this program will require.
*/
//TODO: this also needs to take into account additional input and output maps (once implemented).
-long int estimate_mem_needed ( long int cols, char *mode )
+long int estimate_mem_needed(long int cols, char *mode)
{
- long int mem_count=0;
- long int in_bytes=0;
- long int out_bytes=0;
- long int stat_bytes=0;
+ long int mem_count = 0;
+ long int in_bytes = 0;
+ long int out_bytes = 0;
+ long int stat_bytes = 0;
- /* memory for neighborhood weights and statistics */
- if ( !strcmp (mode, "wmean" )) {
- stat_bytes += sizeof (double) * (DATA_WIDTH*DATA_HEIGHT); /* weights matrix */
- }
- stat_bytes += sizeof (double) * (DATA_WIDTH*DATA_HEIGHT); /* max. cell values */
- stat_bytes += sizeof (double) * (DATA_WIDTH*DATA_HEIGHT); /* max. cell weights */
- stat_bytes += sizeof (unsigned long) * (DATA_WIDTH*DATA_HEIGHT); /* max. cell frequencies */
+ /* memory for neighborhood weights and statistics */
+ if (!strcmp(mode, "wmean")) {
+ stat_bytes += sizeof(double) * (DATA_WIDTH * DATA_HEIGHT); /* weights matrix */
+ }
+ stat_bytes += sizeof(double) * (DATA_WIDTH * DATA_HEIGHT); /* max. cell values */
+ stat_bytes += sizeof(double) * (DATA_WIDTH * DATA_HEIGHT); /* max. cell weights */
+ stat_bytes += sizeof(unsigned long) * (DATA_WIDTH * DATA_HEIGHT); /* max. cell frequencies */
- /* input data rows with padded buffers */
- in_bytes = (unsigned long) (WINDOW_HEIGHT * (cols+(PADDING_WIDTH*2)) );
- if ( IN_TYPE == CELL_TYPE ) {
- in_bytes += in_bytes * sizeof (CELL);
+ /* input data rows with padded buffers */
+ in_bytes = (unsigned long)(WINDOW_HEIGHT * (cols + (PADDING_WIDTH * 2)));
+ if (IN_TYPE == CELL_TYPE) {
+ in_bytes += in_bytes * sizeof(CELL);
}
- if ( IN_TYPE == FCELL_TYPE ) {
- in_bytes += in_bytes * sizeof (FCELL);
+ if (IN_TYPE == FCELL_TYPE) {
+ in_bytes += in_bytes * sizeof(FCELL);
}
- if ( IN_TYPE == DCELL_TYPE ) {
- in_bytes += in_bytes * sizeof (DCELL);
+ if (IN_TYPE == DCELL_TYPE) {
+ in_bytes += in_bytes * sizeof(DCELL);
}
/* output data row */
- out_bytes = (unsigned long) cols;
- if ( OUT_TYPE == CELL_TYPE ) {
- out_bytes += out_bytes * sizeof (CELL);
+ out_bytes = (unsigned long)cols;
+ if (OUT_TYPE == CELL_TYPE) {
+ out_bytes += out_bytes * sizeof(CELL);
}
- if ( OUT_TYPE == FCELL_TYPE ) {
- out_bytes += out_bytes * sizeof (FCELL);
+ if (OUT_TYPE == FCELL_TYPE) {
+ out_bytes += out_bytes * sizeof(FCELL);
}
- if ( OUT_TYPE == DCELL_TYPE ) {
- out_bytes += out_bytes * sizeof (DCELL);
+ if (OUT_TYPE == DCELL_TYPE) {
+ out_bytes += out_bytes * sizeof(DCELL);
}
mem_count = stat_bytes + in_bytes + out_bytes;
- return (mem_count);
+ return (mem_count);
}
@@ -140,35 +142,44 @@
* This uses a fixed layout which may not be able to print very
* large matrices.
*/
-void print_weights_matrix ( long int rows, long int cols )
+void print_weights_matrix(long int rows, long int cols)
{
- int i,j;
- int weight_matrix_line_length = 80;
- int weight_matrix_weight_length = 7;
- char weight_matrix_line_buf[weight_matrix_line_length+1];
- char weight_matrix_weight_buf[weight_matrix_line_length+1];
+ int i, j;
+ int weight_matrix_line_length = 80;
+ int weight_matrix_weight_length = 7;
+ char weight_matrix_line_buf[weight_matrix_line_length + 1];
+ char weight_matrix_weight_buf[weight_matrix_line_length + 1];
- G_message (_("Spatial weights neighborhood (cells):"));
- for ( i = 0; i < rows; i ++ ) {
- G_strncpy (weight_matrix_line_buf, "", 1);
- for ( j = 0; j < cols; j ++ ) {
- if ( WEIGHTS[i][j] != -1.0 ) {
- snprintf (weight_matrix_weight_buf, weight_matrix_line_length, "%06.2f ", WEIGHTS[i][j]);
- if ( strlen (weight_matrix_weight_buf) > (weight_matrix_weight_length) ) {
- snprintf (weight_matrix_weight_buf, weight_matrix_line_length, "[????] ", WEIGHTS[i][j]);
- }
- } else {
- snprintf (weight_matrix_weight_buf, weight_matrix_line_length, "...... ");
- }
- if ( strlen (weight_matrix_weight_buf) + strlen (weight_matrix_line_buf) > weight_matrix_line_length ) {
- G_strncpy ( weight_matrix_line_buf, "[line too long to print]", weight_matrix_line_length);
- break;
- } else {
- G_strcat (weight_matrix_line_buf, weight_matrix_weight_buf);
- }
- }
- fprintf (stdout, "%s\n", weight_matrix_line_buf);
- }
+ G_message(_("Spatial weights neighborhood (cells):"));
+ for (i = 0; i < rows; i++) {
+ G_strncpy(weight_matrix_line_buf, "", 1);
+ for (j = 0; j < cols; j++) {
+ if (WEIGHTS[i][j] != -1.0) {
+ snprintf(weight_matrix_weight_buf, weight_matrix_line_length,
+ "%06.2f ", WEIGHTS[i][j]);
+ if (strlen(weight_matrix_weight_buf) >
+ (weight_matrix_weight_length)) {
+ snprintf(weight_matrix_weight_buf,
+ weight_matrix_line_length, "[????] ",
+ WEIGHTS[i][j]);
+ }
+ }
+ else {
+ snprintf(weight_matrix_weight_buf, weight_matrix_line_length,
+ "...... ");
+ }
+ if (strlen(weight_matrix_weight_buf) +
+ strlen(weight_matrix_line_buf) > weight_matrix_line_length) {
+ G_strncpy(weight_matrix_line_buf, "[line too long to print]",
+ weight_matrix_line_length);
+ break;
+ }
+ else {
+ G_strcat(weight_matrix_line_buf, weight_matrix_weight_buf);
+ }
+ }
+ fprintf(stdout, "%s\n", weight_matrix_line_buf);
+ }
}
@@ -176,19 +187,19 @@
* Returns a void* that points to the first valid data cell in the
* input data row corresponding to "row_idx".
*/
-void *get_input_row ( unsigned long row_idx )
+void *get_input_row(unsigned long row_idx)
{
- unsigned long i;
- void *my_cell=NULL;
+ unsigned long i;
+ void *my_cell = NULL;
- my_cell = CELL_INPUT_HANDLES[row_idx];
+ my_cell = CELL_INPUT_HANDLES[row_idx];
- for ( i=0; i < PADDING_WIDTH; i ++ )
- my_cell += CELL_IN_SIZE;
+ for (i = 0; i < PADDING_WIDTH; i++)
+ my_cell += CELL_IN_SIZE;
- return ( my_cell );
+ return (my_cell);
}
@@ -217,78 +228,88 @@
* types of information to be collected.
*/
-void collect_values_unfiltered ( double val1, double val2, double min, double max, stats_struct *stats )
+void collect_values_unfiltered(double val1, double val2, double min,
+ double max, stats_struct * stats)
{
- stats->values[stats->num_values]=val1;
- stats->certainty+=val2;
- stats->num_values++;
+ stats->values[stats->num_values] = val1;
+ stats->certainty += val2;
+ stats->num_values++;
}
-void collect_values_filtered ( double val1, double val2, double min, double max, stats_struct *stats)
+void collect_values_filtered(double val1, double val2, double min, double max,
+ stats_struct * stats)
{
- unsigned long i;
+ unsigned long i;
- if ( val1 >= min && val1 <= max ) {
- collect_values_unfiltered (val1,val2,min,max,stats);
- }
+ if (val1 >= min && val1 <= max) {
+ collect_values_unfiltered(val1, val2, min, max, stats);
+ }
}
-void collect_values_and_weights_unfiltered ( double val1, double val2, double min, double max, stats_struct *stats )
+void collect_values_and_weights_unfiltered(double val1, double val2,
+ double min, double max,
+ stats_struct * stats)
{
- stats->values[stats->num_values]=val1;
- stats->weights[stats->num_values]=val2;
- stats->certainty+=val2;
- stats->num_values++;
+ stats->values[stats->num_values] = val1;
+ stats->weights[stats->num_values] = val2;
+ stats->certainty += val2;
+ stats->num_values++;
}
-void collect_values_and_weights_filtered ( double val1, double val2, double min, double max, stats_struct *stats )
+void collect_values_and_weights_filtered(double val1, double val2, double min,
+ double max, stats_struct * stats)
{
- unsigned long i;
+ unsigned long i;
- if ( val1 >= min && val1 <= max ) {
- collect_values_and_weights_unfiltered (val1,val2,min,max,stats);
- }
+ if (val1 >= min && val1 <= max) {
+ collect_values_and_weights_unfiltered(val1, val2, min, max, stats);
+ }
}
-void collect_values_and_frequencies_unfiltered ( double val1, double val2, double min, double max, stats_struct *stats )
+void collect_values_and_frequencies_unfiltered(double val1, double val2,
+ double min, double max,
+ stats_struct * stats)
{
- unsigned long i;
- int done;
+ unsigned long i;
+ int done;
- stats->certainty+=val2;
+ stats->certainty += val2;
- /* extreme case: no values collected yet */
- if ( stats->num_values == 0 ) {
- stats->values[0]=val1;
- stats->frequencies[0]=1;
- stats->num_values++;
- return;
- }
+ /* extreme case: no values collected yet */
+ if (stats->num_values == 0) {
+ stats->values[0] = val1;
+ stats->frequencies[0] = 1;
+ stats->num_values++;
+ return;
+ }
- /* at least one value already collected */
- for ( i=0; i < stats->num_values ; i++ ) {
- /* does this value already exist in the stats object? */
- if ( stats->values[i]==val1 ) {
- /* yes: increase its counter and abort search */
- stats->frequencies[i]++;
- stats->values[stats->num_values]=val1;
- stats->num_values++;
- return;
- }
- }
- /* no: first occurrence of this value: store as new entry */
- stats->values[i]=val1;
- stats->frequencies[i]=1;
- stats->num_values++;
+ /* at least one value already collected */
+ for (i = 0; i < stats->num_values; i++) {
+ /* does this value already exist in the stats object? */
+ if (stats->values[i] == val1) {
+ /* yes: increase its counter and abort search */
+ stats->frequencies[i]++;
+ stats->values[stats->num_values] = val1;
+ stats->num_values++;
+ return;
+ }
+ }
+ /* no: first occurrence of this value: store as new entry */
+ stats->values[i] = val1;
+ stats->frequencies[i] = 1;
+ stats->num_values++;
}
-void collect_values_and_frequencies_filtered ( double val1, double val2, double min, double max, stats_struct *stats )
+void collect_values_and_frequencies_filtered(double val1, double val2,
+ double min, double max,
+ stats_struct * stats)
{
- unsigned long i;
+ unsigned long i;
- if ( val1 >= min && val1 <= max ) {
- collect_values_and_frequencies_unfiltered (val1,val2,min,max,stats);
- }
+ if (val1 >= min && val1 <= max) {
+ collect_values_and_frequencies_unfiltered(val1, val2, min, max,
+ stats);
+ }
}
@@ -296,10 +317,13 @@
* Simple double comparision function for use by qsort().
* This is needed for calculating median statistics.
*/
-int compare_dbl (const void* val1, const void* val2) {
- if ( *(double*) val1 == *(double*) val2 ) return 0;
- if ( *(double*) val1 < *(double*) val2 ) return -1;
- return 1;
+int compare_dbl(const void *val1, const void *val2)
+{
+ if (*(double *)val1 == *(double *)val2)
+ return 0;
+ if (*(double *)val1 < *(double *)val2)
+ return -1;
+ return 1;
}
/*
@@ -307,43 +331,48 @@
* basic loop for every type of statistics: collect only non-null
* cells, and collect only within the neighborhood mask.
*/
-void read_neighborhood (unsigned long row_index, unsigned long col,
- double min, double max, int preserve,
- stats_struct *stats )
+void read_neighborhood(unsigned long row_index, unsigned long col,
+ double min, double max, int preserve,
+ stats_struct * stats)
{
- unsigned long i,j;
- void *cell;
- double cell_value;
+ unsigned long i, j;
+ void *cell;
+ double cell_value;
- /* read data */
- unsigned long row_position = row_index-PADDING_HEIGHT;
- stats->num_values=0;
- stats->certainty=0.0;
- for ( i = 0; i < DATA_HEIGHT; i ++ ) {
- cell = CELL_INPUT_HANDLES[i+row_position];
- cell += CELL_IN_SIZE * col;
+ /* read data */
+ unsigned long row_position = row_index - PADDING_HEIGHT;
- for ( j = 0; j < DATA_WIDTH; j ++ ) {
- /* read cell from input buffer */
- if ( !IS_NULL(cell) ) {
- if ( preserve == TRUE && i == (DATA_HEIGHT-1)/2 && j == (DATA_WIDTH-1)/2 ) {
- /* center cell: copy through if needed */
- stats->overwrite = G_malloc (sizeof (double));
- *stats->overwrite = (double)G_get_raster_value_d(cell,IN_TYPE);
- } else {
- /* only add non-null cells to stats */
- cell_value = (double)G_get_raster_value_d(cell,IN_TYPE);
- /* only add if within neighborhood */
- if ( WEIGHTS[i][j] != -1.0 ) {
- /* get data needed for chosen statistic */
- COLLECT_DATA (cell_value, WEIGHTS[i][j],min,max,stats);
- }
- }
- }
- /* go to next cell on current row */
- cell += CELL_IN_SIZE;
- }
- }
+ stats->num_values = 0;
+ stats->certainty = 0.0;
+ for (i = 0; i < DATA_HEIGHT; i++) {
+ cell = CELL_INPUT_HANDLES[i + row_position];
+ cell += CELL_IN_SIZE * col;
+
+ for (j = 0; j < DATA_WIDTH; j++) {
+ /* read cell from input buffer */
+ if (!IS_NULL(cell)) {
+ if (preserve == TRUE && i == (DATA_HEIGHT - 1) / 2 &&
+ j == (DATA_WIDTH - 1) / 2) {
+ /* center cell: copy through if needed */
+ stats->overwrite = G_malloc(sizeof(double));
+ *stats->overwrite =
+ (double)G_get_raster_value_d(cell, IN_TYPE);
+ }
+ else {
+ /* only add non-null cells to stats */
+ cell_value = (double)G_get_raster_value_d(cell, IN_TYPE);
+ /* only add if within neighborhood */
+ if (WEIGHTS[i][j] != -1.0) {
+ /* get data needed for chosen statistic */
+ COLLECT_DATA(cell_value, WEIGHTS[i][j], min, max,
+ stats);
+ }
+ }
+ }
+ /* go to next cell on current row */
+ cell += CELL_IN_SIZE;
+ }
+ }
}
@@ -353,24 +382,24 @@
*
* The cell values are multiplied by their spatial weights before they are stored.
*/
-void get_statistics_wmean (unsigned long row_index, unsigned long col,
- double min, double max, int preserve,
- stats_struct *stats)
+void get_statistics_wmean(unsigned long row_index, unsigned long col,
+ double min, double max, int preserve,
+ stats_struct * stats)
{
- unsigned long i;
- double total;
- double total_weight;
+ unsigned long i;
+ double total;
+ double total_weight;
- read_neighborhood (row_index, col, min, max, preserve, stats );
+ read_neighborhood(row_index, col, min, max, preserve, stats);
- /* compute weighted average of all valid input cells */
- total=0;
- total_weight=0;
- for (i = 0; i < stats->num_values; i ++) {
- total += stats->values[i]*stats->weights[i];
- total_weight += stats->weights[i];
- }
- stats->result = total / total_weight;
+ /* compute weighted average of all valid input cells */
+ total = 0;
+ total_weight = 0;
+ for (i = 0; i < stats->num_values; i++) {
+ total += stats->values[i] * stats->weights[i];
+ total_weight += stats->weights[i];
+ }
+ stats->result = total / total_weight;
}
@@ -378,21 +407,21 @@
* NEIGHBORHOOD STATISTICS FUNCTION MEAN
* Simple, unweighted mean.
*/
-void get_statistics_mean (unsigned long row_index, unsigned long col,
- double min, double max, int preserve,
- stats_struct *stats)
+void get_statistics_mean(unsigned long row_index, unsigned long col,
+ double min, double max, int preserve,
+ stats_struct * stats)
{
- unsigned long i;
- double total;
+ unsigned long i;
+ double total;
- read_neighborhood (row_index, col, min, max, preserve, stats );
+ read_neighborhood(row_index, col, min, max, preserve, stats);
- /* compute total of all valid input cells */
- total=0;
- for (i = 0; i < stats->num_values; i ++) {
- total += stats->values[i];
- }
- stats->result = total / ((double) stats->num_values);
+ /* compute total of all valid input cells */
+ total = 0;
+ for (i = 0; i < stats->num_values; i++) {
+ total += stats->values[i];
+ }
+ stats->result = total / ((double)stats->num_values);
}
@@ -401,22 +430,25 @@
* Simple, unweighted median. For an even number of data points, the median is the
* average of the two central elements in the sorted data list.
*/
-void get_statistics_median (unsigned long row_index, unsigned long col,
- double min, double max, int preserve,
- stats_struct *stats)
+void get_statistics_median(unsigned long row_index, unsigned long col,
+ double min, double max, int preserve,
+ stats_struct * stats)
{
- read_neighborhood (row_index, col, min, max, preserve, stats );
+ read_neighborhood(row_index, col, min, max, preserve, stats);
- /* sort list of values */
- qsort (&stats->values[0], stats->num_values, sizeof (double), &compare_dbl);
+ /* sort list of values */
+ qsort(&stats->values[0], stats->num_values, sizeof(double), &compare_dbl);
- if ( (double)stats->num_values / 2.0 == 0.0 ) {
- /* even number of elements: result is average of the two central values */
- stats->result = (stats->values[stats->num_values/2]+stats->values[(stats->num_values/2)+1])/2.0;
- } else {
- /* odd number of elements: result is the central element */
- stats->result = stats->values[(stats->num_values/2)+1];
- }
+ if ((double)stats->num_values / 2.0 == 0.0) {
+ /* even number of elements: result is average of the two central values */
+ stats->result =
+ (stats->values[stats->num_values / 2] +
+ stats->values[(stats->num_values / 2) + 1]) / 2.0;
+ }
+ else {
+ /* odd number of elements: result is the central element */
+ stats->result = stats->values[(stats->num_values / 2) + 1];
+ }
}
@@ -425,41 +457,41 @@
* Simple, unweighted mode. Mathematically, the mode is not always unique. If there is more than
* one value with highest frequency, the smallest one is chosen to represent the mode.
*/
-void get_statistics_mode (unsigned long row_index, unsigned long col,
- double min, double max, int preserve,
- stats_struct *stats)
+void get_statistics_mode(unsigned long row_index, unsigned long col,
+ double min, double max, int preserve,
+ stats_struct * stats)
{
- unsigned long i;
- double mode;
- unsigned long freq;
+ unsigned long i;
+ double mode;
+ unsigned long freq;
- read_neighborhood ( row_index, col, min, max, preserve, stats );
+ read_neighborhood(row_index, col, min, max, preserve, stats);
- if (stats->num_values < 1)
- return;
+ if (stats->num_values < 1)
+ return;
- mode = stats->values[0];
- freq = stats->frequencies[0];
- for (i=1; i < stats->num_values; i++) {
- if ( stats->frequencies[i] > freq ) {
- mode = stats->values[i];
- freq = stats->frequencies[i];
- }
- }
- stats->result=mode;
+ mode = stats->values[0];
+ freq = stats->frequencies[0];
+ for (i = 1; i < stats->num_values; i++) {
+ if (stats->frequencies[i] > freq) {
+ mode = stats->values[i];
+ freq = stats->frequencies[i];
+ }
+ }
+ stats->result = mode;
}
/*
* Initializes handlers to point to corresponding data rows.
*/
-void init_handles ()
+void init_handles()
{
- unsigned long i;
+ unsigned long i;
- for ( i=0; i < WINDOW_HEIGHT; i ++ ) {
- CELL_INPUT_HANDLES[i] = CELL_INPUT[i];
- }
+ for (i = 0; i < WINDOW_HEIGHT; i++) {
+ CELL_INPUT_HANDLES[i] = CELL_INPUT[i];
+ }
}
@@ -473,88 +505,94 @@
* the data rows for the neighborhood every time we go down one
* row in the current region.
*/
-void advance_one_row (int file_desc, long current_row)
+void advance_one_row(int file_desc, long current_row)
{
- unsigned long i,j;
- void *cell_input;
- static unsigned long replace_row = 0; /* points to the row which will be replaced next */
- unsigned long replace_pos = 0;
+ unsigned long i, j;
+ void *cell_input;
+ static unsigned long replace_row = 0; /* points to the row which will be replaced next */
+ unsigned long replace_pos = 0;
- /* the actual replacement position needs to consider the "no data" padding offset, as well */
- replace_pos = replace_row + PADDING_HEIGHT;
+ /* the actual replacement position needs to consider the "no data" padding offset, as well */
+ replace_pos = replace_row + PADDING_HEIGHT;
- /* get address of data row to replace */
- cell_input = CELL_INPUT[replace_pos];
- for ( i=0; i < PADDING_WIDTH; i ++ )
- cell_input += CELL_IN_SIZE;
+ /* get address of data row to replace */
+ cell_input = CELL_INPUT[replace_pos];
+ for (i = 0; i < PADDING_WIDTH; i++)
+ cell_input += CELL_IN_SIZE;
- /* get next row from disk */
- G_get_raster_row (file_desc, cell_input, current_row+DATA_HEIGHT, IN_TYPE);
+ /* get next row from disk */
+ G_get_raster_row(file_desc, cell_input, current_row + DATA_HEIGHT,
+ IN_TYPE);
- /* re-assign all row handlers below current replacement position */
- j = PADDING_HEIGHT;
- for ( i=0; i < DATA_HEIGHT-(replace_row+1); i ++ ) {
- CELL_INPUT_HANDLES[j] = CELL_INPUT[replace_pos+1+i];
- j++;
- }
+ /* re-assign all row handlers below current replacement position */
+ j = PADDING_HEIGHT;
+ for (i = 0; i < DATA_HEIGHT - (replace_row + 1); i++) {
+ CELL_INPUT_HANDLES[j] = CELL_INPUT[replace_pos + 1 + i];
+ j++;
+ }
- /* re-assign all row handlers up to and including replacement position */
- for ( i=0; i <= replace_row; i ++ ) {
- CELL_INPUT_HANDLES[j] = CELL_INPUT[PADDING_HEIGHT+i];
- j++;
- }
+ /* re-assign all row handlers up to and including replacement position */
+ for (i = 0; i <= replace_row; i++) {
+ CELL_INPUT_HANDLES[j] = CELL_INPUT[PADDING_HEIGHT + i];
+ j++;
+ }
- replace_row ++;
- if ( replace_row > (DATA_HEIGHT-1) ) {
- /* start over once end of data area has been reached */
- replace_row = 0;
- }
+ replace_row++;
+ if (replace_row > (DATA_HEIGHT - 1)) {
+ /* start over once end of data area has been reached */
+ replace_row = 0;
+ }
}
/*
* Interpolates one row of input data, stores result in CELL_OUTPUT (global var)
*/
-void interpolate_row ( unsigned long row_index, unsigned long cols,
- double min, double max, int preserve,
- unsigned long min_cells,
- stats_struct *stats, int write_err )
+void interpolate_row(unsigned long row_index, unsigned long cols,
+ double min, double max, int preserve,
+ unsigned long min_cells,
+ stats_struct * stats, int write_err)
{
- unsigned long j;
- void *cell_input, *cell_output;
- FCELL *err_output;
+ unsigned long j;
+ void *cell_input, *cell_output;
+ FCELL *err_output;
- cell_output = CELL_OUTPUT;
- err_output = ERR_OUTPUT;
+ cell_output = CELL_OUTPUT;
+ err_output = ERR_OUTPUT;
- for ( j = 0; j < cols; j++ ) {
- /* get neighborhood statistics */
- GET_STATS (row_index, j, min, max, preserve, stats);
- /* enough reachable cells in input map? */
- if ( stats->num_values < min_cells ) {
- SET_NULL ( cell_output, 1 );
- if ( write_err )
- G_set_f_null_value (err_output, 1);
- } else {
- if ( stats->overwrite != NULL ) {
- /* write original value into output map */
- WRITE_DOUBLE_VAL (cell_output, *stats->overwrite);
- G_free (stats->overwrite);
- stats->overwrite = NULL;
- } else {
- /* write interpolation result into output map */
- WRITE_DOUBLE_VAL (cell_output, stats->result);
- }
- /* write error/uncertainty output map? */
- if ( write_err ) {
- G_set_raster_value_f(err_output,(FCELL)1.0-(stats->certainty/SUM_WEIGHTS), FCELL_TYPE);
- }
- }
- /* advance cell pointers by one cell size */
- cell_output += CELL_OUT_SIZE;
- err_output ++;
- }
+ for (j = 0; j < cols; j++) {
+ /* get neighborhood statistics */
+ GET_STATS(row_index, j, min, max, preserve, stats);
+ /* enough reachable cells in input map? */
+ if (stats->num_values < min_cells) {
+ SET_NULL(cell_output, 1);
+ if (write_err)
+ G_set_f_null_value(err_output, 1);
+ }
+ else {
+ if (stats->overwrite != NULL) {
+ /* write original value into output map */
+ WRITE_DOUBLE_VAL(cell_output, *stats->overwrite);
+ G_free(stats->overwrite);
+ stats->overwrite = NULL;
+ }
+ else {
+ /* write interpolation result into output map */
+ WRITE_DOUBLE_VAL(cell_output, stats->result);
+ }
+ /* write error/uncertainty output map? */
+ if (write_err) {
+ G_set_raster_value_f(err_output,
+ (FCELL) 1.0 -
+ (stats->certainty / SUM_WEIGHTS),
+ FCELL_TYPE);
+ }
+ }
+ /* advance cell pointers by one cell size */
+ cell_output += CELL_OUT_SIZE;
+ err_output++;
+ }
}
@@ -565,46 +603,50 @@
* passed as "1" and all cells within the circular neighborhood will be set to "1.0".
* In both casses, all cells outside the neighborhood will be set to "-1.0".
*/
-void build_weights_matrix (double radius, double power, double res_x, double res_y, int constant, int use_map_units )
+void build_weights_matrix(double radius, double power, double res_x,
+ double res_y, int constant, int use_map_units)
{
- unsigned long i, j;
- double p1_x, p1_y, p2_x, p2_y, A, B, C, W;
- double tolerance;
+ unsigned long i, j;
+ double p1_x, p1_y, p2_x, p2_y, A, B, C, W;
+ double tolerance;
/* alloc enough mem for weights matrix */
- WEIGHTS=G_malloc (sizeof (double*) * DATA_HEIGHT);
- for ( i = 0; i < DATA_HEIGHT; i ++ ) {
- WEIGHTS[i] = G_malloc (sizeof(double)*DATA_WIDTH);
+ WEIGHTS = G_malloc(sizeof(double *) * DATA_HEIGHT);
+ for (i = 0; i < DATA_HEIGHT; i++) {
+ WEIGHTS[i] = G_malloc(sizeof(double) * DATA_WIDTH);
}
/* center of the neighborhood window in real map units */
- p1_x = (DATA_WIDTH/2 * res_x) + (res_x / 2.0);
- p1_y = (DATA_HEIGHT/2 * res_y) + (res_y / 2.0);
+ p1_x = (DATA_WIDTH / 2 * res_x) + (res_x / 2.0);
+ p1_y = (DATA_HEIGHT / 2 * res_y) + (res_y / 2.0);
/* tolerance for including half cells in the neighborhood */
- tolerance = (sqrt (pow(res_x,2)+pow(res_y,2))) / 2;
+ tolerance = (sqrt(pow(res_x, 2) + pow(res_y, 2))) / 2;
/* 1st pass: get largest possible weight for normalization */
double max = -1.0;
- for ( i = 0; i < DATA_HEIGHT; i ++ ) {
- for ( j = 0; j < DATA_WIDTH; j ++ ) {
+
+ for (i = 0; i < DATA_HEIGHT; i++) {
+ for (j = 0; j < DATA_WIDTH; j++) {
p2_x = (j * res_x) + (res_x / 2.0);
p2_y = (i * res_y) + (res_y / 2.0);
- A = fabs (p2_x - p1_x);
- B = fabs (p2_y - p1_y);
- C = sqrt (pow(A,2)+pow(B,2));
- if ( use_map_units ) {
- if (C > radius+tolerance) {
- WEIGHTS[i][j] = -1.0;
- } else {
- WEIGHTS[i][j] = C;
- }
- } else {
- WEIGHTS[i][j] = C;
- }
- if ( WEIGHTS[i][j] > max ) {
+ A = fabs(p2_x - p1_x);
+ B = fabs(p2_y - p1_y);
+ C = sqrt(pow(A, 2) + pow(B, 2));
+ if (use_map_units) {
+ if (C > radius + tolerance) {
+ WEIGHTS[i][j] = -1.0;
+ }
+ else {
+ WEIGHTS[i][j] = C;
+ }
+ }
+ else {
+ WEIGHTS[i][j] = C;
+ }
+ if (WEIGHTS[i][j] > max) {
/* update max value */
max = WEIGHTS[i][j];
}
@@ -613,43 +655,46 @@
/* Build the weights matrix */
SUM_WEIGHTS = 0.0;
- for ( i = 0; i < DATA_HEIGHT; i ++ ) {
- for ( j = 0; j < DATA_WIDTH; j ++ ) {
+ for (i = 0; i < DATA_HEIGHT; i++) {
+ for (j = 0; j < DATA_WIDTH; j++) {
/* Assign neighborhood coordinates with 0/0
at the top left of the cell neighborhood matrix. */
p2_x = (j * res_x) + (res_x / 2.0);
p2_y = (i * res_y) + (res_y / 2.0);
/* get distance from window center */
- A = fabs (p2_x - p1_x);
- B = fabs (p2_y - p1_y);
- C = sqrt (pow(A,2)+pow(B,2));
- if ( constant ) {
- W = 1.0;
- } else {
- W = ((pow (1-(C/max),power)));
+ A = fabs(p2_x - p1_x);
+ B = fabs(p2_y - p1_y);
+ C = sqrt(pow(A, 2) + pow(B, 2));
+ if (constant) {
+ W = 1.0;
}
+ else {
+ W = ((pow(1 - (C / max), power)));
+ }
/* exclude neighborhood locations that are farther
from the center than the user-defined distance
plus a tolerance of half the current region's
cell diagonal */
- if ( use_map_units ) {
- if (C > radius+tolerance) {
- WEIGHTS[i][j] = -1.0;
- } else {
- WEIGHTS[i][j] = W;
- WEIGHTS[DATA_HEIGHT/2][DATA_WIDTH/2] = 0.0; /* safeguard against total weight growing by center cell weight (InF) */
- SUM_WEIGHTS += WEIGHTS[i][j];
- }
- } else {
- WEIGHTS[i][j] = W;
- WEIGHTS[DATA_HEIGHT/2][DATA_WIDTH/2] = 0.0; /* safeguard against total weight growing by center cell weight (InF) */
- SUM_WEIGHTS += WEIGHTS[i][j];
- }
+ if (use_map_units) {
+ if (C > radius + tolerance) {
+ WEIGHTS[i][j] = -1.0;
+ }
+ else {
+ WEIGHTS[i][j] = W;
+ WEIGHTS[DATA_HEIGHT / 2][DATA_WIDTH / 2] = 0.0; /* safeguard against total weight growing by center cell weight (InF) */
+ SUM_WEIGHTS += WEIGHTS[i][j];
+ }
+ }
+ else {
+ WEIGHTS[i][j] = W;
+ WEIGHTS[DATA_HEIGHT / 2][DATA_WIDTH / 2] = 0.0; /* safeguard against total weight growing by center cell weight (InF) */
+ SUM_WEIGHTS += WEIGHTS[i][j];
+ }
}
}
/* weight of center cell is always = 1 */
- WEIGHTS[DATA_HEIGHT/2][DATA_WIDTH/2] = 1.0;
+ WEIGHTS[DATA_HEIGHT / 2][DATA_WIDTH / 2] = 1.0;
}
@@ -660,8 +705,8 @@
*/
int main(int argc, char *argv[])
{
- /* processing time keepers */
- time_t start, finish;
+ /* processing time keepers */
+ time_t start, finish;
/* GRASS region properties */
struct Cell_head cellhd, region;
@@ -677,11 +722,12 @@
struct GModule *module;
struct
{
- struct Option
- *input, *output, *error,
- *radius, *mode, *power, *min, *max, *minpts;
+ struct Option
+ *input, *output, *error,
+ *radius, *mode, *power, *min, *max, *minpts;
struct Flag
- *dist_m, *preserve, *print_w, *print_u, *center, *single_precision;
+ *dist_m, *preserve, *print_w, *print_u, *center,
+ *single_precision;
} parm;
/* program settings */
@@ -689,13 +735,13 @@
char *output;
char *mapset;
char *error;
- double radius=1.0;
- unsigned long min_cells=12;
- double power=2.0;
- double min=0.0;
- double max=0.0;
- int filter_min=0;
- int filter_max=0;
+ double radius = 1.0;
+ unsigned long min_cells = 12;
+ double power = 2.0;
+ double min = 0.0;
+ double max = 0.0;
+ int filter_min = 0;
+ int filter_max = 0;
int write_error;
/* file handlers */
@@ -705,11 +751,11 @@
int out_fd;
int err_fd;
- /* cell statistics object */
+ /* cell statistics object */
stats_struct cell_stats;
/* generic indices, loop counters, etc. */
- unsigned long i,j,k;
+ unsigned long i, j, k;
long l;
unsigned long count;
@@ -721,7 +767,7 @@
module = G_define_module();
module->keywords = _("raster,interpolation,IDW");
module->description =
- _("Rapidly fills 'no data' cells of a raster map with interpolated values (IDW).");
+ _("Rapidly fills 'no data' cells of a raster map with interpolated values (IDW).");
/* parameters */
@@ -747,7 +793,8 @@
parm.radius->required = YES;
parm.radius->multiple = NO;
parm.radius->type = TYPE_DOUBLE;
- parm.radius->description = _("Distance threshold (default: in cells) for interpolation");
+ parm.radius->description =
+ _("Distance threshold (default: in cells) for interpolation");
parm.radius->answer = "3";
parm.mode = G_define_option();
@@ -766,7 +813,8 @@
parm.min->required = NO;
parm.min->multiple = NO;
parm.min->type = TYPE_DOUBLE;
- parm.min->description = _("Minimum input data value to include in interpolation");
+ parm.min->description =
+ _("Minimum input data value to include in interpolation");
parm.max = G_define_option();
parm.max->key = "maximum";
@@ -774,7 +822,8 @@
parm.max->required = NO;
parm.max->multiple = NO;
parm.max->type = TYPE_DOUBLE;
- parm.max->description = _("Maximum input data value to include in interpolation");
+ parm.max->description =
+ _("Maximum input data value to include in interpolation");
parm.power = G_define_option();
parm.power->key = "power";
@@ -792,11 +841,13 @@
parm.minpts->multiple = NO;
parm.minpts->type = TYPE_INTEGER;
parm.minpts->answer = "8";
- parm.minpts->description = _("Minimum number of data cells within search radius");
+ parm.minpts->description =
+ _("Minimum number of data cells within search radius");
parm.dist_m = G_define_flag();
parm.dist_m->key = 'm';
- parm.dist_m->description = _("Interpret distance as map units, not cell number");
+ parm.dist_m->description =
+ _("Interpret distance as map units, not cell number");
parm.preserve = G_define_flag();
parm.preserve->key = 'p';
@@ -812,11 +863,12 @@
parm.single_precision = G_define_flag();
parm.single_precision->key = 's';
- parm.single_precision->description = _("Single precision floating point output");
+ parm.single_precision->description =
+ _("Single precision floating point output");
if (G_parser(argc, argv))
- exit(EXIT_FAILURE);
+ exit(EXIT_FAILURE);
input = parm.input->answer;
output = parm.output->answer;
@@ -824,211 +876,239 @@
/* get setting of current GRASS region */
G_get_window(®ion);
projection = region.proj;
- if ( projection == PROJECTION_LL && parm.dist_m->answer ) {
- G_warning (_("You are working with lat/lon data."));
- G_warning (_("This module uses a straight-line distance metric."));
- G_warning (_("Expect inaccuracies."));
+ if (projection == PROJECTION_LL && parm.dist_m->answer) {
+ G_warning(_("You are working with lat/lon data."));
+ G_warning(_("This module uses a straight-line distance metric."));
+ G_warning(_("Expect inaccuracies."));
}
- rows = (unsigned long) region.rows;
- cols = (unsigned long) region.cols;
- res_x = (double) region.ew_res;
- res_y = (double) region.ns_res;
+ rows = (unsigned long)region.rows;
+ cols = (unsigned long)region.cols;
+ res_x = (double)region.ew_res;
+ res_y = (double)region.ns_res;
/* get user parameters */
- radius = strtod (parm.radius->answer,0);
- power = strtod (parm.power->answer,0);
- min_cells = atol (parm.minpts->answer);
+ radius = strtod(parm.radius->answer, 0);
+ power = strtod(parm.power->answer, 0);
+ min_cells = atol(parm.minpts->answer);
write_error = 0;
- if ( parm.error->answer ) {
- write_error = 1;
- if (!strcmp(parm.error->answer,parm.output->answer)) {
- G_fatal_error (_("Result map name cannot be identical with uncertainty map name."));
- }
+ if (parm.error->answer) {
+ write_error = 1;
+ if (!strcmp(parm.error->answer, parm.output->answer)) {
+ G_fatal_error(_("Result map name cannot be identical with uncertainty map name."));
+ }
}
/* validate user parameters */
double max_dist = 0.0;
- if ( parm.dist_m->answer ) {
- if ( radius < 0 ) {
- G_fatal_error (_("Maximum distance must be larger than zero."));
+
+ if (parm.dist_m->answer) {
+ if (radius < 0) {
+ G_fatal_error(_("Maximum distance must be larger than zero."));
}
- if ( res_x < res_y ) {
- if ( radius < res_x )
- G_fatal_error (_("Maximum distance must be at least '%.6f' (W-E resolution)."), res_x);
- }
- if ( res_y < res_x ) {
- if ( radius < res_y )
- G_fatal_error (_("Maximum distance must be at least '%.6f' (S-N resolution)."), res_y);
- }
- if ( res_y == res_x ) {
- if ( radius < res_y )
- G_fatal_error (_("Maximum distance must be at least '%.6f' (W-E and S-N resolution)."), res_y);
- }
- max_dist = sqrt ( pow((cols*res_x),2) + pow((rows*res_y),2) );
- if ( radius > max_dist ) {
- G_warning (_("Maximum distance too large. Adjusted to '%.6f' (diagonal of current region)."), max_dist);
- radius = max_dist;
- }
- } else {
+ if (res_x < res_y) {
+ if (radius < res_x)
+ G_fatal_error(_("Maximum distance must be at least '%.6f' (W-E resolution)."),
+ res_x);
+ }
+ if (res_y < res_x) {
+ if (radius < res_y)
+ G_fatal_error(_("Maximum distance must be at least '%.6f' (S-N resolution)."),
+ res_y);
+ }
+ if (res_y == res_x) {
+ if (radius < res_y)
+ G_fatal_error(_("Maximum distance must be at least '%.6f' (W-E and S-N resolution)."),
+ res_y);
+ }
+ max_dist = sqrt(pow((cols * res_x), 2) + pow((rows * res_y), 2));
+ if (radius > max_dist) {
+ G_warning(_("Maximum distance too large. Adjusted to '%.6f' (diagonal of current region)."),
+ max_dist);
+ radius = max_dist;
+ }
+ }
+ else {
unsigned long radius_i = (unsigned long)radius;
- radius=(double)radius_i; //truncate to whole cell number
- if ( radius < 1 ) {
- G_fatal_error (_("Maximum distance must be at least one cell."));
+
+ radius = (double)radius_i; //truncate to whole cell number
+ if (radius < 1) {
+ G_fatal_error(_("Maximum distance must be at least one cell."));
}
- unsigned long max_dist_i = (unsigned long)(sqrt ( pow((cols),2) + pow((rows),2) ));
- max_dist = (double) max_dist_i;
- if ( radius > max_dist ) {
- G_warning (_("Maximum distance too large. Adjusted to '%lu' cells (diagonal of current region)."), max_dist_i);
- radius = (double)max_dist_i;
+ unsigned long max_dist_i =
+ (unsigned long)(sqrt(pow((cols), 2) + pow((rows), 2)));
+ max_dist = (double)max_dist_i;
+ if (radius > max_dist) {
+ G_warning(_("Maximum distance too large. Adjusted to '%lu' cells (diagonal of current region)."),
+ max_dist_i);
+ radius = (double)max_dist_i;
}
}
- if ( min_cells < 1 ) {
- G_fatal_error (_("Minimum number of cells must be at least '1'."));
+ if (min_cells < 1) {
+ G_fatal_error(_("Minimum number of cells must be at least '1'."));
}
- if ( min_cells > ((DATA_WIDTH*DATA_HEIGHT)-1)) {
- G_fatal_error (_("Specified minimum number of cells unreachable with current settings."));
+ if (min_cells > ((DATA_WIDTH * DATA_HEIGHT) - 1)) {
+ G_fatal_error(_("Specified minimum number of cells unreachable with current settings."));
}
- if ( filter_min != 0 && filter_max != 0 ) {
- if ( min >= max ) {
- G_fatal_error (_("Value for 'minimum' must be smaller than value for 'maximum'."));
- }
+ if (filter_min != 0 && filter_max != 0) {
+ if (min >= max) {
+ G_fatal_error(_("Value for 'minimum' must be smaller than value for 'maximum'."));
+ }
}
- if ( parm.power->answer && strcmp (parm.mode->answer,"wmean") ) {
- G_warning (_("The 'power' option has no effect in any mode other than 'wmean'."));
- parm.power->answer=0;
+ if (parm.power->answer && strcmp(parm.mode->answer, "wmean")) {
+ G_warning(_("The 'power' option has no effect in any mode other than 'wmean'."));
+ parm.power->answer = 0;
}
/* rounded dimensions of weight matrix (in map cells) */
- if ( parm.dist_m->answer ) {
- DATA_WIDTH = ((unsigned long) ceil (radius/res_x)) * 2 + 1;
- DATA_HEIGHT = ((unsigned long) ceil (radius/res_y)) * 2 + 1;
- if ((!parm.print_w->answer) && (fmod (radius,res_x) != 0 || fmod (radius,res_y) != 0)) {
- G_warning (_("The specified maximum distance cannot be resolved to whole cells\n at the current resolution settings."));
- }
- } else {
- DATA_WIDTH = (unsigned long) radius*2 + 1;
- DATA_HEIGHT = (unsigned long) radius*2 + 1;
+ if (parm.dist_m->answer) {
+ DATA_WIDTH = ((unsigned long)ceil(radius / res_x)) * 2 + 1;
+ DATA_HEIGHT = ((unsigned long)ceil(radius / res_y)) * 2 + 1;
+ if ((!parm.print_w->answer) &&
+ (fmod(radius, res_x) != 0 || fmod(radius, res_y) != 0)) {
+ G_warning(_("The specified maximum distance cannot be resolved to whole cells\n at the current resolution settings."));
+ }
}
- PADDING_WIDTH = (DATA_WIDTH-1)/2;
- PADDING_HEIGHT = (DATA_HEIGHT-1)/2;
- WINDOW_WIDTH = (PADDING_WIDTH*2) + DATA_WIDTH;
- WINDOW_HEIGHT = (PADDING_HEIGHT*2) + DATA_HEIGHT;
- G_message (_("W-E size of neighborhood is %lu cells."), DATA_WIDTH);
- G_message (_("S-N size of neighborhood is %lu cells."), DATA_HEIGHT);
+ else {
+ DATA_WIDTH = (unsigned long)radius *2 + 1;
+ DATA_HEIGHT = (unsigned long)radius *2 + 1;
+ }
+ PADDING_WIDTH = (DATA_WIDTH - 1) / 2;
+ PADDING_HEIGHT = (DATA_HEIGHT - 1) / 2;
+ WINDOW_WIDTH = (PADDING_WIDTH * 2) + DATA_WIDTH;
+ WINDOW_HEIGHT = (PADDING_HEIGHT * 2) + DATA_HEIGHT;
+ G_message(_("W-E size of neighborhood is %lu cells."), DATA_WIDTH);
+ G_message(_("S-N size of neighborhood is %lu cells."), DATA_HEIGHT);
- if ( DATA_WIDTH < 3 || DATA_HEIGHT < 3 ) {
- G_fatal_error (_("Neighborhood cannot be smaller than 3 cells in X or Y direction."));
+ if (DATA_WIDTH < 3 || DATA_HEIGHT < 3) {
+ G_fatal_error(_("Neighborhood cannot be smaller than 3 cells in X or Y direction."));
}
if (parm.print_w->answer) {
- if ( !strcmp (parm.mode->answer,"wmean") ) {
- build_weights_matrix (radius, power, res_x, res_y, 0, parm.dist_m->answer);
+ if (!strcmp(parm.mode->answer, "wmean")) {
+ build_weights_matrix(radius, power, res_x, res_y, 0,
+ parm.dist_m->answer);
}
- if ( !strcmp (parm.mode->answer,"mean") ) {
- build_weights_matrix (radius, power, res_x, res_y, 1, parm.dist_m->answer);
+ if (!strcmp(parm.mode->answer, "mean")) {
+ build_weights_matrix(radius, power, res_x, res_y, 1,
+ parm.dist_m->answer);
}
- if ( !strcmp (parm.mode->answer,"median") ) {
- build_weights_matrix (radius, power, res_x, res_y, 1, parm.dist_m->answer);
+ if (!strcmp(parm.mode->answer, "median")) {
+ build_weights_matrix(radius, power, res_x, res_y, 1,
+ parm.dist_m->answer);
}
- if ( !strcmp (parm.mode->answer,"mode") ) {
- build_weights_matrix (radius, power, res_x, res_y, 1, parm.dist_m->answer);
+ if (!strcmp(parm.mode->answer, "mode")) {
+ build_weights_matrix(radius, power, res_x, res_y, 1,
+ parm.dist_m->answer);
}
- print_weights_matrix (DATA_HEIGHT, DATA_WIDTH);
+ print_weights_matrix(DATA_HEIGHT, DATA_WIDTH);
/* continue only if "-u" flag has also been given */
- if ( !parm.print_u->answer )
- exit (0);
+ if (!parm.print_u->answer)
+ exit(0);
}
/* open raster input map and get its storage type */
mapset = G_find_cell2(input, "");
if (!mapset)
- G_fatal_error(_("Raster map <%s> not found"), input);
+ G_fatal_error(_("Raster map <%s> not found"), input);
if (G_get_cellhd(input, mapset, &cellhd) < 0)
- G_fatal_error(_("Unable to read header of raster map <%s@%s>"), input, mapset);
+ G_fatal_error(_("Unable to read header of raster map <%s@%s>"), input,
+ mapset);
if ((in_fd = G_open_cell_old(input, mapset)) < 0)
- G_fatal_error(_("Unable to open raster map <%s> in mapset <%s>"), input, mapset);
+ G_fatal_error(_("Unable to open raster map <%s> in mapset <%s>"),
+ input, mapset);
IN_TYPE = G_get_raster_map_type(in_fd);
/* minimum and maximum values for interpolating range */
- if ( IN_TYPE == CELL_TYPE ) {
- G_read_range ( input, mapset, &int_range);
- min = (double) int_range.min;
- max = (double) int_range.max;
- } else {
- G_read_fp_range ( input, mapset, &fp_range);
- min = (double) fp_range.min;
- max = (double) fp_range.max;
+ if (IN_TYPE == CELL_TYPE) {
+ G_read_range(input, mapset, &int_range);
+ min = (double)int_range.min;
+ max = (double)int_range.max;
}
- if ( parm.min->answer ) {
- min = strtod (parm.min->answer,0);
- filter_min = 1;
+ else {
+ G_read_fp_range(input, mapset, &fp_range);
+ min = (double)fp_range.min;
+ max = (double)fp_range.max;
}
- if ( parm.max->answer ) {
- max = strtod (parm.max->answer,0);
- filter_max = 1;
+ if (parm.min->answer) {
+ min = strtod(parm.min->answer, 0);
+ filter_min = 1;
}
- G_message (_("Input data range is %f to %f.\n"), min, max);
+ if (parm.max->answer) {
+ max = strtod(parm.max->answer, 0);
+ filter_max = 1;
+ }
+ G_message(_("Input data range is %f to %f.\n"), min, max);
/* determine input and output data types, advise user */
OUT_TYPE = IN_TYPE;
- if ( IN_TYPE == DCELL_TYPE ) {
- if ( parm.single_precision->answer ) {
- OUT_TYPE = FCELL_TYPE;
- } else {
- OUT_TYPE = DCELL_TYPE;
- }
+ if (IN_TYPE == DCELL_TYPE) {
+ if (parm.single_precision->answer) {
+ OUT_TYPE = FCELL_TYPE;
+ }
+ else {
+ OUT_TYPE = DCELL_TYPE;
+ }
}
- if ( IN_TYPE == CELL_TYPE ) {
- if ( (!strcmp (parm.mode->answer,"wmean")) || (!strcmp (parm.mode->answer,"mean"))
- || (!strcmp (parm.mode->answer,"median"))) {
- G_warning (_("Input data type is integer but interpolation mode is '%s'."), parm.mode->answer);
- if ( parm.single_precision->answer ) {
- OUT_TYPE = FCELL_TYPE;
- G_warning (_("Output type changed to floating point (single)."));
- } else {
- OUT_TYPE = DCELL_TYPE;
- G_warning (_("Output type changed to floating point (double)."));
- }
- } else {
- if ( parm.single_precision->answer ) {
- G_warning (_("Ignoring '%c' flag. Output data type will be integer."), parm.single_precision->key);
- }
- }
+ if (IN_TYPE == CELL_TYPE) {
+ if ((!strcmp(parm.mode->answer, "wmean")) ||
+ (!strcmp(parm.mode->answer, "mean"))
+ || (!strcmp(parm.mode->answer, "median"))) {
+ G_warning(_("Input data type is integer but interpolation mode is '%s'."),
+ parm.mode->answer);
+ if (parm.single_precision->answer) {
+ OUT_TYPE = FCELL_TYPE;
+ G_warning(_("Output type changed to floating point (single)."));
+ }
+ else {
+ OUT_TYPE = DCELL_TYPE;
+ G_warning(_("Output type changed to floating point (double)."));
+ }
+ }
+ else {
+ if (parm.single_precision->answer) {
+ G_warning(_("Ignoring '%c' flag. Output data type will be integer."),
+ parm.single_precision->key);
+ }
+ }
}
char data_type_string_in[21];
char data_type_string_out[21];
- if ( IN_TYPE == CELL_TYPE ) {
- G_strncpy (data_type_string_in, "integer", 20);
- }
- else if ( IN_TYPE == FCELL_TYPE ) {
- G_strncpy (data_type_string_in, "single", 20);
- }
- else if ( IN_TYPE == DCELL_TYPE ) {
- G_strncpy (data_type_string_in, "double", 20);
- }
- if ( OUT_TYPE == CELL_TYPE ) {
- G_strncpy (data_type_string_out, "integer", 20);
- }
- else if ( OUT_TYPE == FCELL_TYPE ) {
- G_strncpy (data_type_string_out, "single", 20);
- }
- else if ( OUT_TYPE == DCELL_TYPE ) {
- G_strncpy (data_type_string_out, "double", 20);
- }
+ if (IN_TYPE == CELL_TYPE) {
+ G_strncpy(data_type_string_in, "integer", 20);
+ }
+ else if (IN_TYPE == FCELL_TYPE) {
+ G_strncpy(data_type_string_in, "single", 20);
+ }
+ else if (IN_TYPE == DCELL_TYPE) {
+ G_strncpy(data_type_string_in, "double", 20);
+ }
+ if (OUT_TYPE == CELL_TYPE) {
+ G_strncpy(data_type_string_out, "integer", 20);
+ }
+ else if (OUT_TYPE == FCELL_TYPE) {
+ G_strncpy(data_type_string_out, "single", 20);
+ }
+ else if (OUT_TYPE == DCELL_TYPE) {
+ G_strncpy(data_type_string_out, "double", 20);
+ }
+
/* initialize data type dependent cell handling functions */
- init_cell_funcs();
+ init_cell_funcs();
- G_message (_("Input data type is '%s' (%i bytes) and output data type is '%s' (%i bytes)."),
- data_type_string_in, CELL_IN_SIZE, data_type_string_out, CELL_OUT_SIZE );
+ G_message(_("Input data type is '%s' (%i bytes) and output data type is '%s' (%i bytes)."),
+ data_type_string_in, CELL_IN_SIZE, data_type_string_out,
+ CELL_OUT_SIZE);
/* just print projected mem usage if user wants it so */
- G_message ("Minimal estimated memory usage is %.3f MB.", ((double) estimate_mem_needed(cols,parm.mode->answer))/1024/1024);
- if ( parm.print_u->answer ) {
- exit (0);
+ G_message("Minimal estimated memory usage is %.3f MB.",
+ ((double)estimate_mem_needed(cols, parm.mode->answer)) / 1024 /
+ 1024);
+ if (parm.print_u->answer) {
+ exit(0);
}
/* Allocate enough memory to read n="distance"x2+1 rows of input map data,
@@ -1038,69 +1118,80 @@
* effect that we can later read data anywhere within the search
* neighborhood, without having to worry about alignment problems.
* */
- CELL_INPUT = G_malloc ( CELL_IN_PTR_SIZE * WINDOW_HEIGHT );
- for ( i = 0; i < WINDOW_HEIGHT; i ++ ) {
- CELL_INPUT[i]= G_malloc ( CELL_IN_SIZE * (cols + (PADDING_WIDTH*2)));
+ CELL_INPUT = G_malloc(CELL_IN_PTR_SIZE * WINDOW_HEIGHT);
+ for (i = 0; i < WINDOW_HEIGHT; i++) {
+ CELL_INPUT[i] = G_malloc(CELL_IN_SIZE * (cols + (PADDING_WIDTH * 2)));
}
- for ( i = 0; i < WINDOW_HEIGHT; i ++ ) {
- G_set_null_value (CELL_INPUT[i], cols + (PADDING_WIDTH*2), IN_TYPE );
+ for (i = 0; i < WINDOW_HEIGHT; i++) {
+ G_set_null_value(CELL_INPUT[i], cols + (PADDING_WIDTH * 2), IN_TYPE);
}
- /*
- * Allocate array of raster row data handlers.
- * When reading rows from the input data buffer, we can use these
- * handlers instead of the original data pointers. That way, we
- * only need to read one new row of data from the disk each time
- * the neighborhood advances down one row in the region. Then, we
- * re-shuffle the row handlers, so that the first handler always
- * points to the first row of the neighborhood data, and the following
- * ones to the subsequent rows.
- * This should save a lot of disk access time.
- */
- CELL_INPUT_HANDLES = G_malloc (CELL_IN_PTR_SIZE*WINDOW_HEIGHT);
+ /*
+ * Allocate array of raster row data handlers.
+ * When reading rows from the input data buffer, we can use these
+ * handlers instead of the original data pointers. That way, we
+ * only need to read one new row of data from the disk each time
+ * the neighborhood advances down one row in the region. Then, we
+ * re-shuffle the row handlers, so that the first handler always
+ * points to the first row of the neighborhood data, and the following
+ * ones to the subsequent rows.
+ * This should save a lot of disk access time.
+ */
+ CELL_INPUT_HANDLES = G_malloc(CELL_IN_PTR_SIZE * WINDOW_HEIGHT);
/* create statistics object */
- cell_stats.values=G_malloc (sizeof(double)*WINDOW_WIDTH*WINDOW_HEIGHT);
- cell_stats.weights=G_malloc (sizeof(double)*WINDOW_WIDTH*WINDOW_HEIGHT);
- cell_stats.frequencies=G_malloc (sizeof(unsigned long)*WINDOW_WIDTH*WINDOW_HEIGHT);
+ cell_stats.values =
+ G_malloc(sizeof(double) * WINDOW_WIDTH * WINDOW_HEIGHT);
+ cell_stats.weights =
+ G_malloc(sizeof(double) * WINDOW_WIDTH * WINDOW_HEIGHT);
+ cell_stats.frequencies =
+ G_malloc(sizeof(unsigned long) * WINDOW_WIDTH * WINDOW_HEIGHT);
cell_stats.overwrite = NULL;
/* set statistics functions according to user option setting */
- if ( !strcmp (parm.mode->answer,"wmean") ) {
- build_weights_matrix (radius, power, res_x, res_y, 0, parm.dist_m->answer);
- GET_STATS = &get_statistics_wmean;
- if ( filter_min == 1 || filter_max == 1 ) {
- COLLECT_DATA=&collect_values_and_weights_filtered;
- } else {
- COLLECT_DATA=&collect_values_and_weights_unfiltered;
- }
+ if (!strcmp(parm.mode->answer, "wmean")) {
+ build_weights_matrix(radius, power, res_x, res_y, 0,
+ parm.dist_m->answer);
+ GET_STATS = &get_statistics_wmean;
+ if (filter_min == 1 || filter_max == 1) {
+ COLLECT_DATA = &collect_values_and_weights_filtered;
+ }
+ else {
+ COLLECT_DATA = &collect_values_and_weights_unfiltered;
+ }
}
- if ( !strcmp (parm.mode->answer,"mean") ) {
- build_weights_matrix (radius, power, res_x, res_y, 1, parm.dist_m->answer);
- GET_STATS = &get_statistics_mean;
- if ( filter_min == 1 || filter_max == 1 ) {
- COLLECT_DATA=&collect_values_filtered;
- } else {
- COLLECT_DATA=&collect_values_unfiltered;
- }
+ if (!strcmp(parm.mode->answer, "mean")) {
+ build_weights_matrix(radius, power, res_x, res_y, 1,
+ parm.dist_m->answer);
+ GET_STATS = &get_statistics_mean;
+ if (filter_min == 1 || filter_max == 1) {
+ COLLECT_DATA = &collect_values_filtered;
+ }
+ else {
+ COLLECT_DATA = &collect_values_unfiltered;
+ }
}
- if ( !strcmp (parm.mode->answer,"median") ) {
- build_weights_matrix (radius, power, res_x, res_y, 1, parm.dist_m->answer);
- GET_STATS = &get_statistics_median;
- if ( filter_min == 1 || filter_max == 1 ) {
- COLLECT_DATA=&collect_values_and_frequencies_filtered;
- } else {
- COLLECT_DATA=&collect_values_and_frequencies_unfiltered;
- }
+ if (!strcmp(parm.mode->answer, "median")) {
+ build_weights_matrix(radius, power, res_x, res_y, 1,
+ parm.dist_m->answer);
+ GET_STATS = &get_statistics_median;
+ if (filter_min == 1 || filter_max == 1) {
+ COLLECT_DATA = &collect_values_and_frequencies_filtered;
+ }
+ else {
+ COLLECT_DATA = &collect_values_and_frequencies_unfiltered;
+ }
}
- if ( !strcmp (parm.mode->answer,"mode") ) {
- build_weights_matrix (radius, power, res_x, res_y, 1, parm.dist_m->answer);
- GET_STATS = &get_statistics_mode;
- if ( filter_min == 1 || filter_max == 1 ) {
- COLLECT_DATA=&collect_values_and_frequencies_filtered;
- } else {
- COLLECT_DATA=&collect_values_and_frequencies_unfiltered;
- }
+ if (!strcmp(parm.mode->answer, "mode")) {
+ build_weights_matrix(radius, power, res_x, res_y, 1,
+ parm.dist_m->answer);
+ GET_STATS = &get_statistics_mode;
+ if (filter_min == 1 || filter_max == 1) {
+ COLLECT_DATA = &collect_values_and_frequencies_filtered;
+ }
+ else {
+ COLLECT_DATA = &collect_values_and_frequencies_unfiltered;
+ }
}
/*
@@ -1112,32 +1203,32 @@
/* Open output map with right data type */
out_fd = G_open_raster_new(output, OUT_TYPE);
if (out_fd < 0) {
- G_fatal_error ("Cannot open output map.");
- exit(EXIT_FAILURE);
+ G_fatal_error("Cannot open output map.");
+ exit(EXIT_FAILURE);
}
/* Reserve memory for one output row buffer */
- CELL_OUTPUT=G_allocate_raster_buf (OUT_TYPE);
+ CELL_OUTPUT = G_allocate_raster_buf(OUT_TYPE);
/* initialize output row */
- SET_NULL (CELL_OUTPUT, cols);
+ SET_NULL(CELL_OUTPUT, cols);
/* produce uncertainty output map? */
- if ( parm.error->answer ) {
- /* Open output map with right data type */
- err_fd = G_open_raster_new(parm.error->answer, FCELL_TYPE);
- if (err_fd < 0) {
- G_fatal_error ("Cannot open uncertainty output map.");
- exit(EXIT_FAILURE);
- }
- ERR_OUTPUT=G_allocate_raster_buf (FCELL_TYPE);
+ if (parm.error->answer) {
+ /* Open output map with right data type */
+ err_fd = G_open_raster_new(parm.error->answer, FCELL_TYPE);
+ if (err_fd < 0) {
+ G_fatal_error("Cannot open uncertainty output map.");
+ exit(EXIT_FAILURE);
+ }
+ ERR_OUTPUT = G_allocate_raster_buf(FCELL_TYPE);
/* initialize output row */
- G_set_f_null_value (ERR_OUTPUT, cols);
+ G_set_f_null_value(ERR_OUTPUT, cols);
}
/* row indices to handle input data buffer */
- unsigned long center_row = (PADDING_HEIGHT*2);
- unsigned long row_idx=0;
+ unsigned long center_row = (PADDING_HEIGHT * 2);
+ unsigned long row_idx = 0;
/* Visit every row in the input dataset.
* To avoid making this code complex and having a lot of
@@ -1147,131 +1238,137 @@
* neighborhood window.
* */
- G_message (_("Interpolating:"));
- unsigned long current_row=0;
+ G_message(_("Interpolating:"));
+ unsigned long current_row = 0;
/* first part: upper edge of region */
- init_handles();
- for ( i=0; i < DATA_HEIGHT; i++ ) {
- cell_input = get_input_row(PADDING_HEIGHT+i);
+ init_handles();
+ for (i = 0; i < DATA_HEIGHT; i++) {
+ cell_input = get_input_row(PADDING_HEIGHT + i);
- cell_input = CELL_INPUT[PADDING_HEIGHT+i];
- for (j=0; j < PADDING_WIDTH; j ++) {
- cell_input += CELL_IN_SIZE;
- }
- G_get_raster_row (in_fd, cell_input, i, IN_TYPE);
- }
- for ( i = 0; i <= PADDING_HEIGHT; i++ ) {
- row_idx = PADDING_HEIGHT + i;
- interpolate_row ( row_idx, cols, min, max, parm.preserve->answer, min_cells,
- &cell_stats, write_error );
- /* write output row buffer to disk */
- G_put_raster_row (out_fd, CELL_OUTPUT, OUT_TYPE);
- if ( parm.error->answer )
- G_put_raster_row (err_fd, ERR_OUTPUT, FCELL_TYPE);
- G_percent(current_row+1, rows, 2);
- current_row++;
+ cell_input = CELL_INPUT[PADDING_HEIGHT + i];
+ for (j = 0; j < PADDING_WIDTH; j++) {
+ cell_input += CELL_IN_SIZE;
+ }
+ G_get_raster_row(in_fd, cell_input, i, IN_TYPE);
}
+ for (i = 0; i <= PADDING_HEIGHT; i++) {
+ row_idx = PADDING_HEIGHT + i;
+ interpolate_row(row_idx, cols, min, max, parm.preserve->answer,
+ min_cells, &cell_stats, write_error);
+ /* write output row buffer to disk */
+ G_put_raster_row(out_fd, CELL_OUTPUT, OUT_TYPE);
+ if (parm.error->answer)
+ G_put_raster_row(err_fd, ERR_OUTPUT, FCELL_TYPE);
+ G_percent(current_row + 1, rows, 2);
+ current_row++;
+ }
/* second part: region between upper and lower edge */
- for ( i = 0; i < rows - (DATA_HEIGHT+1); i++ ) {
- l = i;
- advance_one_row (in_fd,l);
- l++;
- row_idx = center_row;
- interpolate_row ( row_idx, cols, min, max, parm.preserve->answer, min_cells,
- &cell_stats, write_error );
- /* write output row buffer to disk */
- G_put_raster_row (out_fd, CELL_OUTPUT, OUT_TYPE);
- if ( parm.error->answer )
- G_put_raster_row (err_fd, ERR_OUTPUT, FCELL_TYPE);
- G_percent(current_row+1, rows, 2);
- current_row++;
+ for (i = 0; i < rows - (DATA_HEIGHT + 1); i++) {
+ l = i;
+ advance_one_row(in_fd, l);
+ l++;
+ row_idx = center_row;
+ interpolate_row(row_idx, cols, min, max, parm.preserve->answer,
+ min_cells, &cell_stats, write_error);
+ /* write output row buffer to disk */
+ G_put_raster_row(out_fd, CELL_OUTPUT, OUT_TYPE);
+ if (parm.error->answer)
+ G_put_raster_row(err_fd, ERR_OUTPUT, FCELL_TYPE);
+ G_percent(current_row + 1, rows, 2);
+ current_row++;
}
/* third part: lower edge */
init_handles();
- for ( i=rows - DATA_HEIGHT ; i < rows; i++ ) {
- row_idx = (DATA_HEIGHT+PADDING_HEIGHT) - (rows-i);
- cell_input = get_input_row(row_idx);
- G_get_raster_row (in_fd, cell_input, i, IN_TYPE);
- }
- for ( i = rows - PADDING_HEIGHT-1; i < rows; i++ ) {
- row_idx = PADDING_HEIGHT + (DATA_HEIGHT) - (rows-i);
- interpolate_row (row_idx, cols, min, max, parm.preserve->answer, min_cells,
- &cell_stats, write_error );
- /* write output row buffer to disk */
- G_put_raster_row (out_fd, CELL_OUTPUT, OUT_TYPE);
- if ( parm.error->answer )
- G_put_raster_row (err_fd, ERR_OUTPUT, FCELL_TYPE);
- G_percent(current_row+1, rows, 2);
- current_row++;
+ for (i = rows - DATA_HEIGHT; i < rows; i++) {
+ row_idx = (DATA_HEIGHT + PADDING_HEIGHT) - (rows - i);
+ cell_input = get_input_row(row_idx);
+ G_get_raster_row(in_fd, cell_input, i, IN_TYPE);
}
+ for (i = rows - PADDING_HEIGHT - 1; i < rows; i++) {
+ row_idx = PADDING_HEIGHT + (DATA_HEIGHT) - (rows - i);
+ interpolate_row(row_idx, cols, min, max, parm.preserve->answer,
+ min_cells, &cell_stats, write_error);
+ /* write output row buffer to disk */
+ G_put_raster_row(out_fd, CELL_OUTPUT, OUT_TYPE);
+ if (parm.error->answer)
+ G_put_raster_row(err_fd, ERR_OUTPUT, FCELL_TYPE);
+ G_percent(current_row + 1, rows, 2);
+ current_row++;
+ }
/* close all maps */
G_close_cell(out_fd);
G_close_cell(in_fd);
- if ( parm.error->answer ) {
- G_close_cell(err_fd);
+ if (parm.error->answer) {
+ G_close_cell(err_fd);
}
/* Free memory */
- for ( i = 0; i < DATA_HEIGHT; i++ ) {
- G_free (WEIGHTS[i]);
+ for (i = 0; i < DATA_HEIGHT; i++) {
+ G_free(WEIGHTS[i]);
}
- G_free (WEIGHTS);
+ G_free(WEIGHTS);
- if ( CELL_INPUT != NULL ) {
- for ( i = 0; i < DATA_HEIGHT; i ++ ) {
- G_free (CELL_INPUT[i]);
- }
- G_free (CELL_INPUT);
+ if (CELL_INPUT != NULL) {
+ for (i = 0; i < DATA_HEIGHT; i++) {
+ G_free(CELL_INPUT[i]);
+ }
+ G_free(CELL_INPUT);
}
- if ( CELL_OUTPUT != NULL )
- G_free (CELL_OUTPUT);
+ if (CELL_OUTPUT != NULL)
+ G_free(CELL_OUTPUT);
- if ( parm.error->answer )
- G_free (ERR_OUTPUT);
+ if (parm.error->answer)
+ G_free(ERR_OUTPUT);
- G_free (cell_stats.values);
- G_free (cell_stats.weights);
- G_free (cell_stats.frequencies);
+ G_free(cell_stats.values);
+ G_free(cell_stats.weights);
+ G_free(cell_stats.frequencies);
/* write metadata into result and error maps */
- G_short_history (parm.output->answer,"raster",&hist);
- sprintf (hist.title,"Result of interpolation/gap filling");
- hist.edlinecnt=2;
- if ( parm.dist_m->answer ) {
- sprintf (hist.edhist[0],"Settings: mode=%s, distance (map units)=%.6f, power=%.3f",
- parm.mode->answer, radius, power );
- } else {
- sprintf (hist.edhist[0],"Settings: mode=%s, distance (cells)=%lu, power=%.3f",
- parm.mode->answer, (unsigned long) radius, power );
+ G_short_history(parm.output->answer, "raster", &hist);
+ sprintf(hist.title, "Result of interpolation/gap filling");
+ hist.edlinecnt = 2;
+ if (parm.dist_m->answer) {
+ sprintf(hist.edhist[0],
+ "Settings: mode=%s, distance (map units)=%.6f, power=%.3f",
+ parm.mode->answer, radius, power);
}
- sprintf (hist.edhist[1]," min=%.3f, max=%.3f, min. points=%lu",
- min, max, min_cells );
- G_write_history (parm.output->answer, &hist);
+ else {
+ sprintf(hist.edhist[0],
+ "Settings: mode=%s, distance (cells)=%lu, power=%.3f",
+ parm.mode->answer, (unsigned long)radius, power);
+ }
+ sprintf(hist.edhist[1], " min=%.3f, max=%.3f, min. points=%lu",
+ min, max, min_cells);
+ G_write_history(parm.output->answer, &hist);
- if ( parm.error->answer ) {
- G_short_history (parm.error->answer,"raster",&hist);
- sprintf (hist.title,"Uncertainty of interpolation/gap filling");
- hist.edlinecnt=2;
- sprintf (hist.edhist[0],"Result map: %s", parm.output->answer );
- sprintf (hist.edhist[1],"Theoretic range is '0' (lowest) to '1' (highest).");
- G_write_history (parm.error->answer, &hist);
+ if (parm.error->answer) {
+ G_short_history(parm.error->answer, "raster", &hist);
+ sprintf(hist.title, "Uncertainty of interpolation/gap filling");
+ hist.edlinecnt = 2;
+ sprintf(hist.edhist[0], "Result map: %s", parm.output->answer);
+ sprintf(hist.edhist[1],
+ "Theoretic range is '0' (lowest) to '1' (highest).");
+ G_write_history(parm.error->answer, &hist);
}
finish = time(NULL);
- double ticks = difftime(finish,start);
- int hours = trunc (ticks/3600);
- ticks = ticks - (hours*3600);
- int minutes = trunc (ticks/60);
- ticks = ticks - (minutes*60);
- int seconds = trunc (ticks);
+ double ticks = difftime(finish, start);
+ int hours = trunc(ticks / 3600);
+ ticks = ticks - (hours * 3600);
+ int minutes = trunc(ticks / 60);
+
+ ticks = ticks - (minutes * 60);
+ int seconds = trunc(ticks);
+
/* DONE */
- G_done_msg(_("Processing time was %ih%im%is."),hours,minutes,seconds);
+ G_done_msg(_("Processing time was %ih%im%is."), hours, minutes, seconds);
return (EXIT_SUCCESS);
}
More information about the grass-commit
mailing list