[GRASS-SVN] r66595 - grass/trunk/raster/r.in.lidar
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Oct 24 20:42:35 PDT 2015
Author: wenzeslaus
Date: 2015-10-24 20:42:35 -0700 (Sat, 24 Oct 2015)
New Revision: 66595
Added:
grass/trunk/raster/r.in.lidar/point_binning.c
grass/trunk/raster/r.in.lidar/point_binning.h
Modified:
grass/trunk/raster/r.in.lidar/main.c
Log:
r.in.lidar: separate binning code, split to functions
Modified: grass/trunk/raster/r.in.lidar/main.c
===================================================================
--- grass/trunk/raster/r.in.lidar/main.c 2015-10-25 03:10:35 UTC (rev 66594)
+++ grass/trunk/raster/r.in.lidar/main.c 2015-10-25 03:42:35 UTC (rev 66595)
@@ -32,82 +32,18 @@
#include "local_proto.h"
#include "rast_segment.h"
+#include "point_binning.h"
-struct node
-{
- int next;
- double z;
-};
-
-#define SIZE_INCREMENT 10
-int num_nodes = 0;
-int max_nodes = 0;
-struct node *nodes;
-
#define LAS_ALL 0
#define LAS_FIRST 1
#define LAS_LAST 2
#define LAS_MID 3
-int new_node(void)
-{
- int n = num_nodes++;
-
- if (num_nodes >= max_nodes) {
- max_nodes += SIZE_INCREMENT;
- nodes = G_realloc(nodes, (size_t)max_nodes * sizeof(struct node));
- }
-
- return n;
-}
-
-
-/* add node to sorted, single linked list
- * returns id if head has to be saved to index array, otherwise -1 */
-int add_node(int head, double z)
-{
- int node_id, last_id, newnode_id, head_id;
-
- head_id = head;
- node_id = head_id;
- last_id = head_id;
-
- while (node_id != -1 && nodes[node_id].z < z) {
- last_id = node_id;
- node_id = nodes[node_id].next;
- }
-
- /* end of list, simply append */
- if (node_id == -1) {
- newnode_id = new_node();
- nodes[newnode_id].next = -1;
- nodes[newnode_id].z = z;
- nodes[last_id].next = newnode_id;
- return -1;
- }
- else if (node_id == head_id) { /* pole position, insert as head */
- newnode_id = new_node();
- nodes[newnode_id].next = head_id;
- head_id = newnode_id;
- nodes[newnode_id].z = z;
- return (head_id);
- }
- else { /* somewhere in the middle, insert */
- newnode_id = new_node();
- nodes[newnode_id].z = z;
- nodes[newnode_id].next = node_id;
- nodes[last_id].next = newnode_id;
- return -1;
- }
-}
-
int main(int argc, char *argv[])
{
int out_fd, base_raster;
char *infile, *outmap;
int percent;
- int method = -1;
- int bin_n, bin_min, bin_max, bin_sum, bin_sumsq, bin_index;
double zrange_min, zrange_max, d_tmp;
unsigned long estimated_lines;
@@ -115,13 +51,13 @@
struct History history;
char title[64];
SEGMENT base_segment;
- void *n_array, *min_array, *max_array, *sum_array, *sumsq_array,
- *index_array, *base_array;
- void *raster_row, *ptr;
+ struct PointBinning point_binning;
+ void *base_array;
+ void *raster_row;
struct Cell_head region;
struct Cell_head input_region;
int rows, cols; /* scan box size */
- int row, col; /* counters */
+ int row; /* counters */
int pass, npasses;
unsigned long line, line_total;
@@ -134,21 +70,13 @@
int skipme, i;
int point_class;
- double min = 0.0 / 0.0; /* init as nan */
- double max = 0.0 / 0.0; /* init as nan */
double zscale = 1.0;
- size_t offset, n_offset;
- int n = 0;
- double sum = 0.;
- double sumsq = 0.;
- double variance, mean, skew, sumdev;
- int pth = 0;
- double trim = 0.0;
double res = 0.0;
- int j, k;
- int head_id, node_id;
- int r_low, r_up;
+ struct BinIndex bin_index_nodes;
+ bin_index_nodes.num_nodes = 0;
+ bin_index_nodes.max_nodes = 0;
+ bin_index_nodes.nodes = 0;
struct GModule *module;
struct Option *input_opt, *output_opt, *percent_opt, *type_opt, *filter_opt, *class_opt;
@@ -396,104 +324,9 @@
}
}
- /* figure out what maps we need in memory */
- /* n n
- min min
- max max
- range min max max - min
- sum sum
- mean sum n sum/n
- stddev sum sumsq n sqrt((sumsq - sum*sum/n)/n)
- variance sum sumsq n (sumsq - sum*sum/n)/n
- coeff_var sum sumsq n sqrt((sumsq - sum*sum/n)/n) / (sum/n)
- median n array index to linked list
- percentile n array index to linked list
- skewness n array index to linked list
- trimmean n array index to linked list
- */
- bin_n = FALSE;
- bin_min = FALSE;
- bin_max = FALSE;
- bin_sum = FALSE;
- bin_sumsq = FALSE;
- bin_index = FALSE;
+ point_binning_set(&point_binning, method_opt->answer, pth_opt->answer, trim_opt->answer);
- n_array = NULL;
- min_array = NULL;
- max_array = NULL;
- sum_array = NULL;
- sumsq_array = NULL;
- index_array = NULL;
base_array = NULL;
-
- if (strcmp(method_opt->answer, "n") == 0) {
- method = METHOD_N;
- bin_n = TRUE;
- }
- if (strcmp(method_opt->answer, "min") == 0) {
- method = METHOD_MIN;
- bin_min = TRUE;
- }
- if (strcmp(method_opt->answer, "max") == 0) {
- method = METHOD_MAX;
- bin_max = TRUE;
- }
- if (strcmp(method_opt->answer, "range") == 0) {
- method = METHOD_RANGE;
- bin_min = TRUE;
- bin_max = TRUE;
- }
- if (strcmp(method_opt->answer, "sum") == 0) {
- method = METHOD_SUM;
- bin_sum = TRUE;
- }
- if (strcmp(method_opt->answer, "mean") == 0) {
- method = METHOD_MEAN;
- bin_sum = TRUE;
- bin_n = TRUE;
- }
- if (strcmp(method_opt->answer, "stddev") == 0) {
- method = METHOD_STDDEV;
- bin_sum = TRUE;
- bin_sumsq = TRUE;
- bin_n = TRUE;
- }
- if (strcmp(method_opt->answer, "variance") == 0) {
- method = METHOD_VARIANCE;
- bin_sum = TRUE;
- bin_sumsq = TRUE;
- bin_n = TRUE;
- }
- if (strcmp(method_opt->answer, "coeff_var") == 0) {
- method = METHOD_COEFF_VAR;
- bin_sum = TRUE;
- bin_sumsq = TRUE;
- bin_n = TRUE;
- }
- if (strcmp(method_opt->answer, "median") == 0) {
- method = METHOD_MEDIAN;
- bin_index = TRUE;
- }
- if (strcmp(method_opt->answer, "percentile") == 0) {
- if (pth_opt->answer != NULL)
- pth = atoi(pth_opt->answer);
- else
- G_fatal_error(_("Unable to calculate percentile without the pth option specified!"));
- method = METHOD_PERCENTILE;
- bin_index = TRUE;
- }
- if (strcmp(method_opt->answer, "skewness") == 0) {
- method = METHOD_SKEWNESS;
- bin_index = TRUE;
- }
- if (strcmp(method_opt->answer, "trimmean") == 0) {
- if (trim_opt->answer != NULL)
- trim = atof(trim_opt->answer) / 100.0;
- else
- G_fatal_error(_("Unable to calculate trimmed mean without the trim option specified!"));
- method = METHOD_TRIMMEAN;
- bin_index = TRUE;
- }
if (strcmp("CELL", type_opt->answer) == 0)
rtype = CELL_TYPE;
@@ -502,7 +335,7 @@
else
rtype = FCELL_TYPE;
- if (method == METHOD_N)
+ if (point_binning.method == METHOD_N)
rtype = CELL_TYPE;
@@ -584,49 +417,13 @@
}
if (!scan_flag->answer) {
- /* check if rows * (cols + 1) go into a size_t */
- if (sizeof(size_t) < 8) {
- double dsize = rows * (cols + 1);
-
- if (dsize != (size_t)rows * (cols + 1))
+ if (!check_rows_cols_fit_to_size_t(rows, cols))
G_fatal_error(_("Unable to process the hole map at once. "
- "Please set the %s option to some value lower than 100."),
+ "Please set the '%s' option to some value lower than 100."),
percent_opt->key);
+ point_binning_memory_test(&point_binning, rows, cols, rtype);
}
- /* allocate memory (test for enough before we start) */
- if (bin_n)
- n_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
- if (bin_min)
- min_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
- if (bin_max)
- max_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
- if (bin_sum)
- sum_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
- if (bin_sumsq)
- sumsq_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
- if (bin_index)
- index_array =
- G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
- /* TODO: perhaps none of them needs to be freed */
- /* and then free it again */
- if (bin_n)
- G_free(n_array);
- if (bin_min)
- G_free(min_array);
- if (bin_max)
- G_free(max_array);
- if (bin_sum)
- G_free(sum_array);
- if (bin_sumsq)
- G_free(sumsq_array);
- if (bin_index)
- G_free(index_array);
-
- /** end memory test **/
- }
-
-
/* open output map */
out_fd = Rast_open_new(outmap, rtype);
@@ -672,39 +469,8 @@
G_debug(2, "pass=%d/%d pass_n=%f pass_s=%f rows=%d",
pass, npasses, pass_north, pass_south, rows);
+ point_binning_allocate(&point_binning, rows, cols, rtype);
- if (bin_n) {
- G_debug(2, "allocating n_array");
- n_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
- blank_array(n_array, rows, cols, CELL_TYPE, 0);
- }
- if (bin_min) {
- G_debug(2, "allocating min_array");
- min_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
- blank_array(min_array, rows, cols, rtype, -1); /* fill with NULLs */
- }
- if (bin_max) {
- G_debug(2, "allocating max_array");
- max_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
- blank_array(max_array, rows, cols, rtype, -1); /* fill with NULLs */
- }
- if (bin_sum) {
- G_debug(2, "allocating sum_array");
- sum_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
- blank_array(sum_array, rows, cols, rtype, 0);
- }
- if (bin_sumsq) {
- G_debug(2, "allocating sumsq_array");
- sumsq_array = G_calloc((size_t)rows * (cols + 1), Rast_cell_size(rtype));
- blank_array(sumsq_array, rows, cols, rtype, 0);
- }
- if (bin_index) {
- G_debug(2, "allocating index_array");
- index_array =
- G_calloc((size_t)rows * (cols + 1), Rast_cell_size(CELL_TYPE));
- blank_array(index_array, rows, cols, CELL_TYPE, -1); /* fill with NULLs */
- }
-
line = 0;
count = 0;
counter = 0;
@@ -814,37 +580,7 @@
count++;
/* G_debug(5, "x: %f, y: %f, z: %f", x, y, z); */
- if (bin_n)
- update_n(n_array, cols, arr_row, arr_col);
- if (bin_min)
- update_min(min_array, cols, arr_row, arr_col, rtype, z);
- if (bin_max)
- update_max(max_array, cols, arr_row, arr_col, rtype, z);
- if (bin_sum)
- update_sum(sum_array, cols, arr_row, arr_col, rtype, z);
- if (bin_sumsq)
- update_sumsq(sumsq_array, cols, arr_row, arr_col, rtype, z);
- if (bin_index) {
- ptr = index_array;
- ptr =
- G_incr_void_ptr(ptr,
- ((arr_row * cols) +
- arr_col) * Rast_cell_size(CELL_TYPE));
-
- if (Rast_is_null_value(ptr, CELL_TYPE)) { /* first node */
- head_id = new_node();
- nodes[head_id].next = -1;
- nodes[head_id].z = z;
- Rast_set_c_value(ptr, head_id, CELL_TYPE); /* store index to head */
- }
- else { /* head is already there */
-
- head_id = Rast_get_c_value(ptr, CELL_TYPE); /* get index to head */
- head_id = add_node(head_id, z);
- if (head_id != -1)
- Rast_set_c_value(ptr, head_id, CELL_TYPE); /* store index to head */
- }
- }
+ update_value(&point_binning, &bin_index_nodes, cols, arr_row, arr_col, rtype, z);
} /* while !EOF */
G_percent(1, 1, 1); /* flush */
@@ -855,327 +591,13 @@
/* calc stats and output */
G_message(_("Writing to map ..."));
for (row = 0; row < rows; row++) {
-
- switch (method) {
- case METHOD_N: /* n is a straight copy */
- Rast_raster_cpy(raster_row,
- n_array +
- (row * cols * Rast_cell_size(CELL_TYPE)), cols,
- CELL_TYPE);
- break;
-
- case METHOD_MIN:
- Rast_raster_cpy(raster_row,
- min_array + (row * cols * Rast_cell_size(rtype)),
- cols, rtype);
- break;
-
- case METHOD_MAX:
- Rast_raster_cpy(raster_row,
- max_array + (row * cols * Rast_cell_size(rtype)),
- cols, rtype);
- break;
-
- case METHOD_SUM:
- Rast_raster_cpy(raster_row,
- sum_array + (row * cols * Rast_cell_size(rtype)),
- cols, rtype);
- break;
-
- case METHOD_RANGE: /* (max-min) */
- ptr = raster_row;
- for (col = 0; col < cols; col++) {
- offset = (row * cols + col) * Rast_cell_size(rtype);
- min = Rast_get_d_value(min_array + offset, rtype);
- max = Rast_get_d_value(max_array + offset, rtype);
- Rast_set_d_value(ptr, max - min, rtype);
- ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
- }
- break;
-
- case METHOD_MEAN: /* (sum / n) */
- ptr = raster_row;
- for (col = 0; col < cols; col++) {
- offset = (row * cols + col) * Rast_cell_size(rtype);
- n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
- n = Rast_get_c_value(n_array + n_offset, CELL_TYPE);
- sum = Rast_get_d_value(sum_array + offset, rtype);
-
- if (n == 0)
- Rast_set_null_value(ptr, 1, rtype);
- else
- Rast_set_d_value(ptr, (sum / n), rtype);
-
- ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
- }
- break;
-
- case METHOD_STDDEV: /* sqrt(variance) */
- case METHOD_VARIANCE: /* (sumsq - sum*sum/n)/n */
- case METHOD_COEFF_VAR: /* 100 * stdev / mean */
- ptr = raster_row;
- for (col = 0; col < cols; col++) {
- offset = (row * cols + col) * Rast_cell_size(rtype);
- n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
- n = Rast_get_c_value(n_array + n_offset, CELL_TYPE);
- sum = Rast_get_d_value(sum_array + offset, rtype);
- sumsq = Rast_get_d_value(sumsq_array + offset, rtype);
-
- if (n == 0)
- Rast_set_null_value(ptr, 1, rtype);
- else if (n == 1)
- Rast_set_d_value(ptr, 0.0, rtype);
- else {
- variance = (sumsq - sum * sum / n) / n;
- if (variance < GRASS_EPSILON)
- variance = 0.0;
-
- /* nan test */
- if (variance != variance)
- Rast_set_null_value(ptr, 1, rtype);
- else {
-
- if (method == METHOD_STDDEV)
- variance = sqrt(variance);
-
- else if (method == METHOD_COEFF_VAR)
- variance = 100 * sqrt(variance) / (sum / n);
-
- /* nan test */
- if (variance != variance)
- variance = 0.0; /* OK for n > 0 ?*/
-
- Rast_set_d_value(ptr, variance, rtype);
- }
-
- }
- ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
- }
- break;
- case METHOD_MEDIAN: /* median, if only one point in cell we will use that */
- ptr = raster_row;
- for (col = 0; col < cols; col++) {
- n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
- if (Rast_is_null_value(index_array + n_offset, CELL_TYPE)) /* no points in cell */
- Rast_set_null_value(ptr, 1, rtype);
- else { /* one or more points in cell */
-
- head_id =
- Rast_get_c_value(index_array + n_offset,
- CELL_TYPE);
- node_id = head_id;
-
- n = 0;
-
- while (node_id != -1) { /* count number of points in cell */
- n++;
- node_id = nodes[node_id].next;
- }
-
- if (n == 1) /* only one point, use that */
- Rast_set_d_value(ptr, nodes[head_id].z,
- rtype);
- else if (n % 2 != 0) { /* odd number of points: median_i = (n + 1) / 2 */
- n = (n + 1) / 2;
- node_id = head_id;
- for (j = 1; j < n; j++) /* get "median element" */
- node_id = nodes[node_id].next;
-
- Rast_set_d_value(ptr, nodes[node_id].z,
- rtype);
- }
- else { /* even number of points: median = (val_below + val_above) / 2 */
-
- z = (n + 1) / 2.0;
- n = floor(z);
- node_id = head_id;
- for (j = 1; j < n; j++) /* get element "below" */
- node_id = nodes[node_id].next;
-
- z = (nodes[node_id].z +
- nodes[nodes[node_id].next].z) / 2;
- Rast_set_d_value(ptr, z, rtype);
- }
- }
- ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
- }
- break;
- case METHOD_PERCENTILE: /* rank = (pth*(n+1))/100; interpolate linearly */
- ptr = raster_row;
- for (col = 0; col < cols; col++) {
- n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
- if (Rast_is_null_value(index_array + n_offset, CELL_TYPE)) /* no points in cell */
- Rast_set_null_value(ptr, 1, rtype);
- else {
- head_id =
- Rast_get_c_value(index_array + n_offset,
- CELL_TYPE);
- node_id = head_id;
- n = 0;
-
- while (node_id != -1) { /* count number of points in cell */
- n++;
- node_id = nodes[node_id].next;
- }
-
- z = (pth * (n + 1)) / 100.0;
- r_low = floor(z); /* lower rank */
- if (r_low < 1)
- r_low = 1;
- else if (r_low > n)
- r_low = n;
-
- r_up = ceil(z); /* upper rank */
- if (r_up > n)
- r_up = n;
-
- node_id = head_id;
- for (j = 1; j < r_low; j++) /* search lower value */
- node_id = nodes[node_id].next;
-
- z = nodes[node_id].z; /* save lower value */
- node_id = head_id;
- for (j = 1; j < r_up; j++) /* search upper value */
- node_id = nodes[node_id].next;
-
- z = (z + nodes[node_id].z) / 2;
- Rast_set_d_value(ptr, z, rtype);
- }
- ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
- }
- break;
- case METHOD_SKEWNESS: /* skewness = sum(xi-mean)^3/(N-1)*s^3 */
- ptr = raster_row;
- for (col = 0; col < cols; col++) {
- n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
- if (Rast_is_null_value(index_array + n_offset, CELL_TYPE)) /* no points in cell */
- Rast_set_null_value(ptr, 1, rtype);
- else {
- head_id =
- Rast_get_c_value(index_array + n_offset,
- CELL_TYPE);
- node_id = head_id;
-
- n = 0; /* count */
- sum = 0.0; /* sum */
- sumsq = 0.0; /* sum of squares */
- sumdev = 0.0; /* sum of (xi - mean)^3 */
- skew = 0.0; /* skewness */
-
- while (node_id != -1) {
- z = nodes[node_id].z;
- n++;
- sum += z;
- sumsq += (z * z);
- node_id = nodes[node_id].next;
- }
-
- if (n > 1) { /* if n == 1, skew is "0.0" */
- mean = sum / n;
- node_id = head_id;
- while (node_id != -1) {
- z = nodes[node_id].z;
- sumdev += pow((z - mean), 3);
- node_id = nodes[node_id].next;
- }
-
- variance = (sumsq - sum * sum / n) / n;
- if (variance < GRASS_EPSILON)
- skew = 0.0;
- else
- skew =
- sumdev / ((n - 1) *
- pow(sqrt(variance), 3));
- }
- Rast_set_d_value(ptr, skew, rtype);
- }
- ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
- }
- break;
- case METHOD_TRIMMEAN:
- ptr = raster_row;
- for (col = 0; col < cols; col++) {
- n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
- if (Rast_is_null_value(index_array + n_offset, CELL_TYPE)) /* no points in cell */
- Rast_set_null_value(ptr, 1, rtype);
- else {
- head_id =
- Rast_get_c_value(index_array + n_offset,
- CELL_TYPE);
-
- node_id = head_id;
- n = 0;
- while (node_id != -1) { /* count number of points in cell */
- n++;
- node_id = nodes[node_id].next;
- }
-
- if (1 == n)
- mean = nodes[head_id].z;
- else {
- k = floor(trim * n + 0.5); /* number of ranks to discard on each tail */
-
- if (k > 0 && (n - 2 * k) > 0) { /* enough elements to discard */
- node_id = head_id;
- for (j = 0; j < k; j++) /* move to first rank to consider */
- node_id = nodes[node_id].next;
-
- j = k + 1;
- k = n - k;
- n = 0;
- sum = 0.0;
-
- while (j <= k) { /* get values in interval */
- n++;
- sum += nodes[node_id].z;
- node_id = nodes[node_id].next;
- j++;
- }
- }
- else {
- node_id = head_id;
- n = 0;
- sum = 0.0;
- while (node_id != -1) {
- n++;
- sum += nodes[node_id].z;
- node_id = nodes[node_id].next;
- }
- }
- mean = sum / n;
- }
- Rast_set_d_value(ptr, mean, rtype);
- }
- ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
- }
- break;
-
- default:
- G_fatal_error("?");
- }
-
+ write_values(&point_binning, &bin_index_nodes, raster_row, row, cols, rtype);
/* write out line of raster data */
Rast_put_row(out_fd, raster_row, rtype);
}
/* free memory */
- if (bin_n)
- G_free(n_array);
- if (bin_min)
- G_free(min_array);
- if (bin_max)
- G_free(max_array);
- if (bin_sum)
- G_free(sum_array);
- if (bin_sumsq)
- G_free(sumsq_array);
- if (bin_index) {
- G_free(index_array);
- G_free(nodes);
- num_nodes = 0;
- max_nodes = 0;
- nodes = NULL;
- }
+ point_binning_free(&point_binning, &bin_index_nodes);
} /* passes loop */
if (base_array)
Rast_close(base_raster);
Added: grass/trunk/raster/r.in.lidar/point_binning.c
===================================================================
--- grass/trunk/raster/r.in.lidar/point_binning.c (rev 0)
+++ grass/trunk/raster/r.in.lidar/point_binning.c 2015-10-25 03:42:35 UTC (rev 66595)
@@ -0,0 +1,735 @@
+/*
+ * r.in.lidar projection-related functions
+ *
+ * Copyright 2011-2015 by Markus Metz, and The GRASS Development Team
+ * Authors:
+ * Markus Metz (r.in.lidar)
+ * Vaclav Petras (move code to separate functions)
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/raster.h>
+
+#include "point_binning.h"
+#include "local_proto.h"
+
+#define SIZE_INCREMENT 10
+
+static int new_node(struct BinIndex *bin_index)
+{
+ int n = bin_index->num_nodes++;
+
+ if (bin_index->num_nodes >= bin_index->max_nodes) {
+ bin_index->max_nodes += SIZE_INCREMENT;
+ bin_index->nodes = G_realloc(bin_index->nodes,
+ (size_t) bin_index->max_nodes *
+ sizeof(struct node));
+ }
+
+ return n;
+}
+
+
+/* add node to sorted, single linked list
+ * returns id if head has to be saved to index array, otherwise -1 */
+static int add_node(struct BinIndex *bin_index, int head, double z)
+{
+ int node_id, last_id, newnode_id, head_id;
+
+ head_id = head;
+ node_id = head_id;
+ last_id = head_id;
+
+ while (node_id != -1 && bin_index->nodes[node_id].z < z) {
+ last_id = node_id;
+ node_id = bin_index->nodes[node_id].next;
+ }
+
+ /* end of list, simply append */
+ if (node_id == -1) {
+ newnode_id = new_node(bin_index);
+ bin_index->nodes[newnode_id].next = -1;
+ bin_index->nodes[newnode_id].z = z;
+ bin_index->nodes[last_id].next = newnode_id;
+ return -1;
+ }
+ else if (node_id == head_id) { /* pole position, insert as head */
+ newnode_id = new_node(bin_index);
+ bin_index->nodes[newnode_id].next = head_id;
+ head_id = newnode_id;
+ bin_index->nodes[newnode_id].z = z;
+ return (head_id);
+ }
+ else { /* somewhere in the middle, insert */
+ newnode_id = new_node(bin_index);
+ bin_index->nodes[newnode_id].z = z;
+ bin_index->nodes[newnode_id].next = node_id;
+ bin_index->nodes[last_id].next = newnode_id;
+ return -1;
+ }
+}
+
+
+int update_bin_index(struct BinIndex *bin_index, void *index_array,
+ int cols, int row, int col,
+ RASTER_MAP_TYPE map_type, double value)
+{
+ int head_id;
+ void *ptr = index_array;
+
+ ptr =
+ G_incr_void_ptr(ptr,
+ ((row * cols) + col) * Rast_cell_size(CELL_TYPE));
+
+ /* first node */
+ if (Rast_is_null_value(ptr, CELL_TYPE)) {
+ head_id = new_node(bin_index);
+ bin_index->nodes[head_id].next = -1;
+ bin_index->nodes[head_id].z = value;
+ /* store index to head */
+ Rast_set_c_value(ptr, head_id, CELL_TYPE);
+ }
+ /* head is already there */
+ else {
+
+ /* get index to head */
+ head_id = Rast_get_c_value(ptr, CELL_TYPE);
+ head_id = add_node(bin_index, head_id, value);
+ /* if id valid, store index to head */
+ if (head_id != -1)
+ Rast_set_c_value(ptr, head_id, CELL_TYPE);
+ }
+ /* for consistency with functions from support.c */
+ return 0;
+}
+
+void point_binning_set(struct PointBinning *point_binning, char *method,
+ char *percentile, char *trim)
+{
+
+ /* figure out what maps we need in memory */
+ /* n n
+ min min
+ max max
+ range min max max - min
+ sum sum
+ mean sum n sum/n
+ stddev sum sumsq n sqrt((sumsq - sum*sum/n)/n)
+ variance sum sumsq n (sumsq - sum*sum/n)/n
+ coeff_var sum sumsq n sqrt((sumsq - sum*sum/n)/n) / (sum/n)
+ median n array index to linked list
+ percentile n array index to linked list
+ skewness n array index to linked list
+ trimmean n array index to linked list
+ */
+ point_binning->bin_n = FALSE;
+ point_binning->bin_min = FALSE;
+ point_binning->bin_max = FALSE;
+ point_binning->bin_sum = FALSE;
+ point_binning->bin_sumsq = FALSE;
+ point_binning->bin_index = FALSE;
+
+ point_binning->n_array = NULL;
+ point_binning->min_array = NULL;
+ point_binning->max_array = NULL;
+ point_binning->sum_array = NULL;
+ point_binning->sumsq_array = NULL;
+ point_binning->index_array = NULL;
+
+ if (strcmp(method, "n") == 0) {
+ point_binning->method = METHOD_N;
+ point_binning->bin_n = TRUE;
+ }
+ if (strcmp(method, "min") == 0) {
+ point_binning->method = METHOD_MIN;
+ point_binning->bin_min = TRUE;
+ }
+ if (strcmp(method, "max") == 0) {
+ point_binning->method = METHOD_MAX;
+ point_binning->bin_max = TRUE;
+ }
+ if (strcmp(method, "range") == 0) {
+ point_binning->method = METHOD_RANGE;
+ point_binning->bin_min = TRUE;
+ point_binning->bin_max = TRUE;
+ }
+ if (strcmp(method, "sum") == 0) {
+ point_binning->method = METHOD_SUM;
+ point_binning->bin_sum = TRUE;
+ }
+ if (strcmp(method, "mean") == 0) {
+ point_binning->method = METHOD_MEAN;
+ point_binning->bin_sum = TRUE;
+ point_binning->bin_n = TRUE;
+ }
+ if (strcmp(method, "stddev") == 0) {
+ point_binning->method = METHOD_STDDEV;
+ point_binning->bin_sum = TRUE;
+ point_binning->bin_sumsq = TRUE;
+ point_binning->bin_n = TRUE;
+ }
+ if (strcmp(method, "variance") == 0) {
+ point_binning->method = METHOD_VARIANCE;
+ point_binning->bin_sum = TRUE;
+ point_binning->bin_sumsq = TRUE;
+ point_binning->bin_n = TRUE;
+ }
+ if (strcmp(method, "coeff_var") == 0) {
+ point_binning->method = METHOD_COEFF_VAR;
+ point_binning->bin_sum = TRUE;
+ point_binning->bin_sumsq = TRUE;
+ point_binning->bin_n = TRUE;
+ }
+ if (strcmp(method, "median") == 0) {
+ point_binning->method = METHOD_MEDIAN;
+ point_binning->bin_index = TRUE;
+ }
+ if (strcmp(method, "percentile") == 0) {
+ if (percentile != NULL)
+ point_binning->pth = atoi(percentile);
+ else
+ G_fatal_error(_("Unable to calculate percentile without the pth option specified!"));
+ point_binning->method = METHOD_PERCENTILE;
+ point_binning->bin_index = TRUE;
+ }
+ if (strcmp(method, "skewness") == 0) {
+ point_binning->method = METHOD_SKEWNESS;
+ point_binning->bin_index = TRUE;
+ }
+ if (strcmp(method, "trimmean") == 0) {
+ if (trim != NULL)
+ point_binning->trim = atof(trim) / 100.0;
+ else
+ G_fatal_error(_("Unable to calculate trimmed mean without the trim option specified!"));
+ point_binning->method = METHOD_TRIMMEAN;
+ point_binning->bin_index = TRUE;
+ }
+}
+
+
+/* check if rows * (cols + 1) go into a size_t */
+int check_rows_cols_fit_to_size_t(int rows, int cols)
+{
+ if (sizeof(size_t) < 8) {
+ double dsize = rows * (cols + 1);
+
+ if (dsize != (size_t) rows * (cols + 1))
+ return FALSE;
+ }
+ return TRUE;
+}
+
+
+void point_binning_memory_test(struct PointBinning *point_binning, int rows,
+ int cols, RASTER_MAP_TYPE rtype)
+{
+ /* allocate memory (test for enough before we start) */
+ if (point_binning->bin_n)
+ point_binning->n_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(CELL_TYPE));
+ if (point_binning->bin_min)
+ point_binning->min_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
+ if (point_binning->bin_max)
+ point_binning->max_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
+ if (point_binning->bin_sum)
+ point_binning->sum_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
+ if (point_binning->bin_sumsq)
+ point_binning->sumsq_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
+ if (point_binning->bin_index)
+ point_binning->index_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(CELL_TYPE));
+ /* TODO: perhaps none of them needs to be freed */
+
+ /* and then free it again */
+ if (point_binning->bin_n)
+ G_free(point_binning->n_array);
+ if (point_binning->bin_min)
+ G_free(point_binning->min_array);
+ if (point_binning->bin_max)
+ G_free(point_binning->max_array);
+ if (point_binning->bin_sum)
+ G_free(point_binning->sum_array);
+ if (point_binning->bin_sumsq)
+ G_free(point_binning->sumsq_array);
+ if (point_binning->bin_index)
+ G_free(point_binning->index_array);
+}
+
+
+void point_binning_allocate(struct PointBinning *point_binning, int rows,
+ int cols, RASTER_MAP_TYPE rtype)
+{
+ if (point_binning->bin_n) {
+ G_debug(2, "allocating n_array");
+ point_binning->n_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(CELL_TYPE));
+ blank_array(point_binning->n_array, rows, cols, CELL_TYPE, 0);
+ }
+ if (point_binning->bin_min) {
+ G_debug(2, "allocating min_array");
+ point_binning->min_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
+ blank_array(point_binning->min_array, rows, cols, rtype, -1); /* fill with NULLs */
+ }
+ if (point_binning->bin_max) {
+ G_debug(2, "allocating max_array");
+ point_binning->max_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
+ blank_array(point_binning->max_array, rows, cols, rtype, -1); /* fill with NULLs */
+ }
+ if (point_binning->bin_sum) {
+ G_debug(2, "allocating sum_array");
+ point_binning->sum_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
+ blank_array(point_binning->sum_array, rows, cols, rtype, 0);
+ }
+ if (point_binning->bin_sumsq) {
+ G_debug(2, "allocating sumsq_array");
+ point_binning->sumsq_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
+ blank_array(point_binning->sumsq_array, rows, cols, rtype, 0);
+ }
+ if (point_binning->bin_index) {
+ G_debug(2, "allocating index_array");
+ point_binning->index_array =
+ G_calloc((size_t) rows * (cols + 1), Rast_cell_size(CELL_TYPE));
+ blank_array(point_binning->index_array, rows, cols, CELL_TYPE, -1); /* fill with NULLs */
+ }
+}
+
+void point_binning_free(struct PointBinning *point_binning,
+ struct BinIndex *bin_index_nodes)
+{
+ if (point_binning->bin_n)
+ G_free(point_binning->n_array);
+ if (point_binning->bin_min)
+ G_free(point_binning->min_array);
+ if (point_binning->bin_max)
+ G_free(point_binning->max_array);
+ if (point_binning->bin_sum)
+ G_free(point_binning->sum_array);
+ if (point_binning->bin_sumsq)
+ G_free(point_binning->sumsq_array);
+ if (point_binning->bin_index) {
+ G_free(point_binning->index_array);
+ G_free(bin_index_nodes->nodes);
+ bin_index_nodes->num_nodes = 0;
+ bin_index_nodes->max_nodes = 0;
+ bin_index_nodes->nodes = NULL;
+ }
+}
+
+void write_variance(void *raster_row, void *n_array, void *sum_array,
+ void *sumsq_array, int row, int cols,
+ RASTER_MAP_TYPE rtype, int method)
+{
+ size_t offset, n_offset;
+ int n = 0;
+ double variance;
+ double sum = 0.;
+ double sumsq = 0.;
+ int col;
+ void *ptr = raster_row;
+
+ for (col = 0; col < cols; col++) {
+ offset = (row * cols + col) * Rast_cell_size(rtype);
+ n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
+ n = Rast_get_c_value(n_array + n_offset, CELL_TYPE);
+ sum = Rast_get_d_value(sum_array + offset, rtype);
+ sumsq = Rast_get_d_value(sumsq_array + offset, rtype);
+
+ if (n == 0)
+ Rast_set_null_value(ptr, 1, rtype);
+ else if (n == 1)
+ Rast_set_d_value(ptr, 0.0, rtype);
+ else {
+ variance = (sumsq - sum * sum / n) / n;
+ if (variance < GRASS_EPSILON)
+ variance = 0.0;
+
+ /* nan test */
+ if (variance != variance)
+ Rast_set_null_value(ptr, 1, rtype);
+ else {
+
+ if (method == METHOD_STDDEV)
+ variance = sqrt(variance);
+
+ else if (method == METHOD_COEFF_VAR)
+ variance = 100 * sqrt(variance) / (sum / n);
+
+ /* nan test */
+ if (variance != variance)
+ variance = 0.0; /* OK for n > 0 ? */
+
+ Rast_set_d_value(ptr, variance, rtype);
+ }
+
+ }
+ ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+ }
+}
+
+void write_median(struct BinIndex *bin_index, void *raster_row,
+ void *index_array, int row, int cols, RASTER_MAP_TYPE rtype)
+{
+ size_t n_offset;
+ int n;
+ int j;
+ double z;
+ int col;
+ int node_id, head_id;
+ void *ptr = raster_row;
+
+ for (col = 0; col < cols; col++) {
+ n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
+ if (Rast_is_null_value(index_array + n_offset, CELL_TYPE)) /* no points in cell */
+ Rast_set_null_value(ptr, 1, rtype);
+ else { /* one or more points in cell */
+
+ head_id = Rast_get_c_value(index_array + n_offset, CELL_TYPE);
+ node_id = head_id;
+
+ n = 0;
+
+ while (node_id != -1) { /* count number of points in cell */
+ n++;
+ node_id = bin_index->nodes[node_id].next;
+ }
+
+ if (n == 1) /* only one point, use that */
+ Rast_set_d_value(ptr, bin_index->nodes[head_id].z, rtype);
+ else if (n % 2 != 0) { /* odd number of points: median_i = (n + 1) / 2 */
+ n = (n + 1) / 2;
+ node_id = head_id;
+ for (j = 1; j < n; j++) /* get "median element" */
+ node_id = bin_index->nodes[node_id].next;
+
+ Rast_set_d_value(ptr, bin_index->nodes[node_id].z, rtype);
+ }
+ else { /* even number of points: median = (val_below + val_above) / 2 */
+
+ z = (n + 1) / 2.0;
+ n = floor(z);
+ node_id = head_id;
+ for (j = 1; j < n; j++) /* get element "below" */
+ node_id = bin_index->nodes[node_id].next;
+
+ z = (bin_index->nodes[node_id].z +
+ bin_index->nodes[bin_index->nodes[node_id].next].z) / 2;
+ Rast_set_d_value(ptr, z, rtype);
+ }
+ }
+ ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+ }
+}
+
+void write_percentile(struct BinIndex *bin_index, void *raster_row,
+ void *index_array, int row, int cols,
+ RASTER_MAP_TYPE rtype, int pth)
+{
+ size_t n_offset;
+ int n;
+ int j;
+ double z;
+ int col;
+ int node_id, head_id;
+ int r_low, r_up;
+ void *ptr = raster_row;
+
+ for (col = 0; col < cols; col++) {
+ n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
+ if (Rast_is_null_value(index_array + n_offset, CELL_TYPE)) /* no points in cell */
+ Rast_set_null_value(ptr, 1, rtype);
+ else {
+ head_id = Rast_get_c_value(index_array + n_offset, CELL_TYPE);
+ node_id = head_id;
+ n = 0;
+
+ while (node_id != -1) { /* count number of points in cell */
+ n++;
+ node_id = bin_index->nodes[node_id].next;
+ }
+
+ z = (pth * (n + 1)) / 100.0;
+ r_low = floor(z); /* lower rank */
+ if (r_low < 1)
+ r_low = 1;
+ else if (r_low > n)
+ r_low = n;
+
+ r_up = ceil(z); /* upper rank */
+ if (r_up > n)
+ r_up = n;
+
+ node_id = head_id;
+ for (j = 1; j < r_low; j++) /* search lower value */
+ node_id = bin_index->nodes[node_id].next;
+
+ z = bin_index->nodes[node_id].z; /* save lower value */
+ node_id = head_id;
+ for (j = 1; j < r_up; j++) /* search upper value */
+ node_id = bin_index->nodes[node_id].next;
+
+ z = (z + bin_index->nodes[node_id].z) / 2;
+ Rast_set_d_value(ptr, z, rtype);
+ }
+ ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+ }
+}
+
+void write_skewness(struct BinIndex *bin_index, void *raster_row,
+ void *index_array, int row, int cols,
+ RASTER_MAP_TYPE rtype)
+{
+ size_t n_offset;
+ int n;
+ double z;
+ int col;
+ int node_id, head_id;
+ double variance, mean, skew, sumdev;
+ double sum = 0.;
+ double sumsq = 0.;
+ void *ptr = raster_row;
+
+ for (col = 0; col < cols; col++) {
+ n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
+ if (Rast_is_null_value(index_array + n_offset, CELL_TYPE)) /* no points in cell */
+ Rast_set_null_value(ptr, 1, rtype);
+ else {
+ head_id = Rast_get_c_value(index_array + n_offset, CELL_TYPE);
+ node_id = head_id;
+
+ n = 0; /* count */
+ sum = 0.0; /* sum */
+ sumsq = 0.0; /* sum of squares */
+ sumdev = 0.0; /* sum of (xi - mean)^3 */
+ skew = 0.0; /* skewness */
+
+ while (node_id != -1) {
+ z = bin_index->nodes[node_id].z;
+ n++;
+ sum += z;
+ sumsq += (z * z);
+ node_id = bin_index->nodes[node_id].next;
+ }
+
+ if (n > 1) { /* if n == 1, skew is "0.0" */
+ mean = sum / n;
+ node_id = head_id;
+ while (node_id != -1) {
+ z = bin_index->nodes[node_id].z;
+ sumdev += pow((z - mean), 3);
+ node_id = bin_index->nodes[node_id].next;
+ }
+
+ variance = (sumsq - sum * sum / n) / n;
+ if (variance < GRASS_EPSILON)
+ skew = 0.0;
+ else
+ skew = sumdev / ((n - 1) * pow(sqrt(variance), 3));
+ }
+ Rast_set_d_value(ptr, skew, rtype);
+ }
+ ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+ }
+}
+
+void write_trimmean(struct BinIndex *bin_index, void *raster_row,
+ void *index_array, int row, int cols,
+ RASTER_MAP_TYPE rtype, double trim)
+{
+ size_t n_offset;
+ int n;
+ int j, k;
+ int col;
+ int node_id, head_id;
+ double mean;
+ double sum = 0.;
+ void *ptr = raster_row;
+
+ for (col = 0; col < cols; col++) {
+ n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
+ if (Rast_is_null_value(index_array + n_offset, CELL_TYPE)) /* no points in cell */
+ Rast_set_null_value(ptr, 1, rtype);
+ else {
+ head_id = Rast_get_c_value(index_array + n_offset, CELL_TYPE);
+
+ node_id = head_id;
+ n = 0;
+ while (node_id != -1) { /* count number of points in cell */
+ n++;
+ node_id = bin_index->nodes[node_id].next;
+ }
+
+ if (1 == n)
+ mean = bin_index->nodes[head_id].z;
+ else {
+ k = floor(trim * n + 0.5); /* number of ranks to discard on each tail */
+
+ if (k > 0 && (n - 2 * k) > 0) { /* enough elements to discard */
+ node_id = head_id;
+ for (j = 0; j < k; j++) /* move to first rank to consider */
+ node_id = bin_index->nodes[node_id].next;
+
+ j = k + 1;
+ k = n - k;
+ n = 0;
+ sum = 0.0;
+
+ while (j <= k) { /* get values in interval */
+ n++;
+ sum += bin_index->nodes[node_id].z;
+ node_id = bin_index->nodes[node_id].next;
+ j++;
+ }
+ }
+ else {
+ node_id = head_id;
+ n = 0;
+ sum = 0.0;
+ while (node_id != -1) {
+ n++;
+ sum += bin_index->nodes[node_id].z;
+ node_id = bin_index->nodes[node_id].next;
+ }
+ }
+ mean = sum / n;
+ }
+ Rast_set_d_value(ptr, mean, rtype);
+ }
+ ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+ }
+}
+
+void write_values(struct PointBinning *point_binning,
+ struct BinIndex *bin_index_nodes, void *raster_row, int row,
+ int cols, RASTER_MAP_TYPE rtype)
+{
+ void *ptr = NULL;
+ int col;
+
+ switch (point_binning->method) {
+ case METHOD_N: /* n is a straight copy */
+ Rast_raster_cpy(raster_row,
+ point_binning->n_array +
+ (row * cols * Rast_cell_size(CELL_TYPE)), cols,
+ CELL_TYPE);
+ break;
+
+ case METHOD_MIN:
+ Rast_raster_cpy(raster_row,
+ point_binning->min_array +
+ (row * cols * Rast_cell_size(rtype)), cols, rtype);
+ break;
+
+ case METHOD_MAX:
+ Rast_raster_cpy(raster_row,
+ point_binning->max_array +
+ (row * cols * Rast_cell_size(rtype)), cols, rtype);
+ break;
+
+ case METHOD_SUM:
+ Rast_raster_cpy(raster_row,
+ point_binning->sum_array +
+ (row * cols * Rast_cell_size(rtype)), cols, rtype);
+ break;
+
+ case METHOD_RANGE: /* (max-min) */
+ ptr = raster_row;
+ for (col = 0; col < cols; col++) {
+ size_t offset = (row * cols + col) * Rast_cell_size(rtype);
+ double min =
+ Rast_get_d_value(point_binning->min_array + offset, rtype);
+ double max =
+ Rast_get_d_value(point_binning->max_array + offset, rtype);
+ Rast_set_d_value(ptr, max - min, rtype);
+ ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+ }
+ break;
+
+ case METHOD_MEAN: /* (sum / n) */
+ ptr = raster_row;
+ for (col = 0; col < cols; col++) {
+ size_t offset = (row * cols + col) * Rast_cell_size(rtype);
+ size_t n_offset = (row * cols + col) * Rast_cell_size(CELL_TYPE);
+ int n = Rast_get_c_value(point_binning->n_array + n_offset,
+ CELL_TYPE);
+ double sum =
+ Rast_get_d_value(point_binning->sum_array + offset, rtype);
+
+ if (n == 0)
+ Rast_set_null_value(ptr, 1, rtype);
+ else
+ Rast_set_d_value(ptr, (sum / n), rtype);
+
+ ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
+ }
+ break;
+
+ case METHOD_STDDEV: /* sqrt(variance) */
+ case METHOD_VARIANCE: /* (sumsq - sum*sum/n)/n */
+ case METHOD_COEFF_VAR: /* 100 * stdev / mean */
+ write_variance(raster_row, point_binning->n_array,
+ point_binning->sum_array, point_binning->sumsq_array,
+ row, cols, rtype, point_binning->method);
+ break;
+ case METHOD_MEDIAN: /* median, if only one point in cell we will use that */
+ write_median(bin_index_nodes, raster_row, point_binning->index_array,
+ row, cols, rtype);
+ break;
+ case METHOD_PERCENTILE: /* rank = (pth*(n+1))/100; interpolate linearly */
+ write_percentile(bin_index_nodes, raster_row,
+ point_binning->index_array, row, cols, rtype,
+ point_binning->pth);
+ break;
+ case METHOD_SKEWNESS: /* skewness = sum(xi-mean)^3/(N-1)*s^3 */
+ write_skewness(bin_index_nodes, raster_row,
+ point_binning->index_array, row, cols, rtype);
+ break;
+ case METHOD_TRIMMEAN:
+ write_trimmean(bin_index_nodes, raster_row,
+ point_binning->index_array, row, cols, rtype,
+ point_binning->trim);
+ break;
+
+ default:
+ G_fatal_error("?");
+ }
+}
+
+void update_value(struct PointBinning *point_binning,
+ struct BinIndex *bin_index_nodes, int cols, int arr_row,
+ int arr_col, RASTER_MAP_TYPE rtype, double z)
+{
+ if (point_binning->bin_n)
+ update_n(point_binning->n_array, cols, arr_row, arr_col);
+ if (point_binning->bin_min)
+ update_min(point_binning->min_array, cols, arr_row, arr_col, rtype,
+ z);
+ if (point_binning->bin_max)
+ update_max(point_binning->max_array, cols, arr_row, arr_col, rtype,
+ z);
+ if (point_binning->bin_sum)
+ update_sum(point_binning->sum_array, cols, arr_row, arr_col, rtype,
+ z);
+ if (point_binning->bin_sumsq)
+ update_sumsq(point_binning->sumsq_array, cols, arr_row, arr_col,
+ rtype, z);
+ if (point_binning->bin_index)
+ update_bin_index(bin_index_nodes, point_binning->index_array, cols,
+ arr_row, arr_col, rtype, z);
+}
Property changes on: grass/trunk/raster/r.in.lidar/point_binning.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass/trunk/raster/r.in.lidar/point_binning.h
===================================================================
--- grass/trunk/raster/r.in.lidar/point_binning.h (rev 0)
+++ grass/trunk/raster/r.in.lidar/point_binning.h 2015-10-25 03:42:35 UTC (rev 66595)
@@ -0,0 +1,95 @@
+/*
+ * r.in.lidar projection-related functions
+ *
+ * Copyright 2011-2015 by Markus Metz, and The GRASS Development Team
+ * Authors:
+ * Markus Metz (r.in.lidar)
+ * Vaclav Petras (move code to separate functions)
+ *
+ * This program is free software licensed under the GPL (>=v2).
+ * Read the COPYING file that comes with GRASS for details.
+ *
+ */
+
+#ifndef __POINT_BINNING_H__
+#define __POINT_BINNING_H__
+
+
+#include <grass/raster.h>
+
+
+struct node
+{
+ int next;
+ double z;
+};
+
+struct BinIndex
+{
+ int num_nodes;
+ int max_nodes;
+ struct node *nodes;
+};
+
+struct PointBinning
+{
+ int method;
+
+ int bin_n;
+ int bin_min;
+ int bin_max;
+ int bin_sum;
+ int bin_sumsq;
+ int bin_index;
+
+ void *n_array;
+ void *min_array;
+ void *max_array;
+ void *sum_array;
+ void *sumsq_array;
+ void *index_array;
+
+ int pth;
+ double trim;
+};
+
+int check_rows_cols_fit_to_size_t(int rows, int cols);
+void point_binning_memory_test(struct PointBinning *point_binning, int rows,
+ int cols, RASTER_MAP_TYPE rtype);
+
+void point_binning_set(struct PointBinning *point_binning, char *method,
+ char *percentile, char *trim);
+void point_binning_allocate(struct PointBinning *point_binning, int rows,
+ int cols, RASTER_MAP_TYPE rtype);
+
+void point_binning_free(struct PointBinning *point_binning,
+ struct BinIndex *bin_index_nodes);
+
+int update_bin_index(struct BinIndex *bin_index, void *index_array,
+ int cols, int row, int col,
+ RASTER_MAP_TYPE map_type, double value);
+void write_variance(void *raster_row, void *n_array, void *sum_array,
+ void *sumsq_array, int row, int cols,
+ RASTER_MAP_TYPE rtype, int method);
+void write_median(struct BinIndex *bin_index, void *raster_row,
+ void *index_array, int row, int cols,
+ RASTER_MAP_TYPE rtype);
+void write_percentile(struct BinIndex *bin_index, void *raster_row,
+ void *index_array, int row, int cols,
+ RASTER_MAP_TYPE rtype, int pth);
+void write_skewness(struct BinIndex *bin_index, void *raster_row,
+ void *index_array, int row, int cols,
+ RASTER_MAP_TYPE rtype);
+void write_trimmean(struct BinIndex *bin_index, void *raster_row,
+ void *index_array, int row, int cols,
+ RASTER_MAP_TYPE rtype, double trim);
+
+void write_values(struct PointBinning *point_binning,
+ struct BinIndex *bin_index_nodes, void *raster_row, int row,
+ int cols, RASTER_MAP_TYPE rtype);
+void update_value(struct PointBinning *point_binning,
+ struct BinIndex *bin_index_nodes, int cols, int arr_row,
+ int arr_col, RASTER_MAP_TYPE rtype, double z);
+
+
+#endif /* __POINT_BINNING_H__ */
Property changes on: grass/trunk/raster/r.in.lidar/point_binning.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list