1 /*
2 * r.in.lidar projection-related functions
3 *
4 * Copyright 2011-2015 by Markus Metz, and The GRASS Development Team
5 * Authors:
6 * Markus Metz (r.in.lidar)
7 * Vaclav Petras (move code to separate functions)
8 *
9 * This program is free software licensed under the GPL (>=v2).
10 * Read the COPYING file that comes with GRASS for details.
11 *
12 */
13
14 #include <stdlib.h>
15 #include <string.h>
16 #include <math.h>
17
18 #include <grass/gis.h>
19 #include <grass/glocale.h>
20 #include <grass/raster.h>
21 #include <grass/vector.h>
22
23 #include "point_binning.h"
24 #include "local_proto.h"
25
26 #define SIZE_INCREMENT 10
27
new_node(struct BinIndex * bin_index)28 static int new_node(struct BinIndex *bin_index)
29 {
30 int n = bin_index->num_nodes++;
31
32 if (bin_index->num_nodes >= bin_index->max_nodes) {
33 bin_index->max_nodes += SIZE_INCREMENT;
34 bin_index->nodes = G_realloc(bin_index->nodes,
35 (size_t) bin_index->max_nodes *
36 sizeof(struct node));
37 }
38
39 return n;
40 }
41
42
43 /* add node to sorted, single linked list
44 * returns id if head has to be saved to index array, otherwise -1 */
add_node(struct BinIndex * bin_index,int head,double z)45 static int add_node(struct BinIndex *bin_index, int head, double z)
46 {
47 int node_id, last_id, newnode_id, head_id;
48
49 head_id = head;
50 node_id = head_id;
51 last_id = head_id;
52
53 while (node_id != -1 && bin_index->nodes[node_id].z < z) {
54 last_id = node_id;
55 node_id = bin_index->nodes[node_id].next;
56 }
57
58 /* end of list, simply append */
59 if (node_id == -1) {
60 newnode_id = new_node(bin_index);
61 bin_index->nodes[newnode_id].next = -1;
62 bin_index->nodes[newnode_id].z = z;
63 bin_index->nodes[last_id].next = newnode_id;
64 return -1;
65 }
66 else if (node_id == head_id) { /* pole position, insert as head */
67 newnode_id = new_node(bin_index);
68 bin_index->nodes[newnode_id].next = head_id;
69 head_id = newnode_id;
70 bin_index->nodes[newnode_id].z = z;
71 return (head_id);
72 }
73 else { /* somewhere in the middle, insert */
74 newnode_id = new_node(bin_index);
75 bin_index->nodes[newnode_id].z = z;
76 bin_index->nodes[newnode_id].next = node_id;
77 bin_index->nodes[last_id].next = newnode_id;
78 return -1;
79 }
80 }
81
82
update_bin_index(struct BinIndex * bin_index,void * index_array,int cols,int row,int col,RASTER_MAP_TYPE map_type,double value)83 int update_bin_index(struct BinIndex *bin_index, void *index_array,
84 int cols, int row, int col,
85 RASTER_MAP_TYPE map_type, double value)
86 {
87 int head_id;
88 void *ptr = index_array;
89
90 ptr =
91 G_incr_void_ptr(ptr,
92 (((size_t) row * cols) + col) * Rast_cell_size(CELL_TYPE));
93
94 /* first node */
95 if (Rast_is_null_value(ptr, CELL_TYPE)) {
96 head_id = new_node(bin_index);
97 bin_index->nodes[head_id].next = -1;
98 bin_index->nodes[head_id].z = value;
99 /* store index to head */
100 Rast_set_c_value(ptr, head_id, CELL_TYPE);
101 }
102 /* head is already there */
103 else {
104
105 /* get index to head */
106 head_id = Rast_get_c_value(ptr, CELL_TYPE);
107 head_id = add_node(bin_index, head_id, value);
108 /* if id valid, store index to head */
109 if (head_id != -1)
110 Rast_set_c_value(ptr, head_id, CELL_TYPE);
111 }
112 /* for consistency with functions from support.c */
113 return 0;
114 }
115
point_binning_set(struct PointBinning * point_binning,char * method,char * percentile,char * trim,int bin_coordinates)116 void point_binning_set(struct PointBinning *point_binning, char *method,
117 char *percentile, char *trim, int bin_coordinates)
118 {
119
120 /* figure out what maps we need in memory */
121 /* n n
122 min min
123 max max
124 range min max max - min
125 sum sum
126 mean sum n sum/n
127 stddev sum sumsq n sqrt((sumsq - sum*sum/n)/n)
128 variance sum sumsq n (sumsq - sum*sum/n)/n
129 coeff_var sum sumsq n sqrt((sumsq - sum*sum/n)/n) / (sum/n)
130 median n array index to linked list
131 percentile n array index to linked list
132 skewness n array index to linked list
133 trimmean n array index to linked list
134 */
135 point_binning->method = METHOD_NONE;
136 point_binning->bin_n = FALSE;
137 point_binning->bin_min = FALSE;
138 point_binning->bin_max = FALSE;
139 point_binning->bin_sum = FALSE;
140 point_binning->bin_sumsq = FALSE;
141 point_binning->bin_index = FALSE;
142 point_binning->bin_coordinates = FALSE;
143
144 point_binning->n_array = NULL;
145 point_binning->min_array = NULL;
146 point_binning->max_array = NULL;
147 point_binning->sum_array = NULL;
148 point_binning->sumsq_array = NULL;
149 point_binning->index_array = NULL;
150 point_binning->x_array = NULL;
151 point_binning->y_array = NULL;
152
153 if (strcmp(method, "n") == 0) {
154 point_binning->method = METHOD_N;
155 point_binning->bin_n = TRUE;
156 }
157 if (strcmp(method, "min") == 0) {
158 point_binning->method = METHOD_MIN;
159 point_binning->bin_min = TRUE;
160 }
161 if (strcmp(method, "max") == 0) {
162 point_binning->method = METHOD_MAX;
163 point_binning->bin_max = TRUE;
164 }
165 if (strcmp(method, "range") == 0) {
166 point_binning->method = METHOD_RANGE;
167 point_binning->bin_min = TRUE;
168 point_binning->bin_max = TRUE;
169 }
170 if (strcmp(method, "sum") == 0) {
171 point_binning->method = METHOD_SUM;
172 point_binning->bin_sum = TRUE;
173 }
174 if (strcmp(method, "mean") == 0) {
175 point_binning->method = METHOD_MEAN;
176 point_binning->bin_sum = TRUE;
177 point_binning->bin_n = TRUE;
178 }
179 if (strcmp(method, "stddev") == 0) {
180 point_binning->method = METHOD_STDDEV;
181 point_binning->bin_sum = TRUE;
182 point_binning->bin_sumsq = TRUE;
183 point_binning->bin_n = TRUE;
184 }
185 if (strcmp(method, "variance") == 0) {
186 point_binning->method = METHOD_VARIANCE;
187 point_binning->bin_sum = TRUE;
188 point_binning->bin_sumsq = TRUE;
189 point_binning->bin_n = TRUE;
190 }
191 if (strcmp(method, "coeff_var") == 0) {
192 point_binning->method = METHOD_COEFF_VAR;
193 point_binning->bin_sum = TRUE;
194 point_binning->bin_sumsq = TRUE;
195 point_binning->bin_n = TRUE;
196 }
197 if (strcmp(method, "median") == 0) {
198 point_binning->method = METHOD_MEDIAN;
199 point_binning->bin_index = TRUE;
200 }
201 if (strcmp(method, "percentile") == 0) {
202 if (percentile != NULL)
203 point_binning->pth = atoi(percentile);
204 else
205 G_fatal_error(_("Unable to calculate percentile without the pth option specified!"));
206 point_binning->method = METHOD_PERCENTILE;
207 point_binning->bin_index = TRUE;
208 }
209 if (strcmp(method, "skewness") == 0) {
210 point_binning->method = METHOD_SKEWNESS;
211 point_binning->bin_index = TRUE;
212 }
213 if (strcmp(method, "trimmean") == 0) {
214 if (trim != NULL)
215 point_binning->trim = atof(trim) / 100.0;
216 else
217 G_fatal_error(_("Unable to calculate trimmed mean without the trim option specified!"));
218 point_binning->method = METHOD_TRIMMEAN;
219 point_binning->bin_index = TRUE;
220 }
221 if (bin_coordinates) {
222 /* x, y */
223 point_binning->bin_coordinates = TRUE;
224 /* z and n */
225 point_binning->bin_sum = TRUE;
226 point_binning->bin_n = TRUE;
227 }
228 }
229
230
231 /* check if rows * (cols + 1) go into a size_t */
check_rows_cols_fit_to_size_t(int rows,int cols)232 int check_rows_cols_fit_to_size_t(int rows, int cols)
233 {
234 if (sizeof(size_t) < 8) {
235 double dsize = rows * (cols + 1);
236
237 /* TODO: the comparison with double may fail */
238 if (dsize != (size_t) rows * (cols + 1))
239 return FALSE;
240 }
241 return TRUE;
242 }
243
244
point_binning_memory_test(struct PointBinning * point_binning,int rows,int cols,RASTER_MAP_TYPE rtype)245 void point_binning_memory_test(struct PointBinning *point_binning, int rows,
246 int cols, RASTER_MAP_TYPE rtype)
247 {
248 /* allocate memory (test for enough before we start) */
249 if (point_binning->bin_n)
250 point_binning->n_array =
251 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(CELL_TYPE));
252 if (point_binning->bin_min)
253 point_binning->min_array =
254 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
255 if (point_binning->bin_max)
256 point_binning->max_array =
257 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
258 if (point_binning->bin_sum)
259 point_binning->sum_array =
260 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
261 if (point_binning->bin_sumsq)
262 point_binning->sumsq_array =
263 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
264 if (point_binning->bin_index)
265 point_binning->index_array =
266 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(CELL_TYPE));
267 if (point_binning->bin_coordinates) {
268 point_binning->x_array =
269 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
270 point_binning->y_array =
271 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
272 }
273 /* TODO: perhaps none of them needs to be freed */
274
275 /* and then free it again */
276 if (point_binning->bin_n)
277 G_free(point_binning->n_array);
278 if (point_binning->bin_min)
279 G_free(point_binning->min_array);
280 if (point_binning->bin_max)
281 G_free(point_binning->max_array);
282 if (point_binning->bin_sum)
283 G_free(point_binning->sum_array);
284 if (point_binning->bin_sumsq)
285 G_free(point_binning->sumsq_array);
286 if (point_binning->bin_index)
287 G_free(point_binning->index_array);
288 if (point_binning->bin_coordinates) {
289 G_free(point_binning->x_array);
290 G_free(point_binning->y_array);
291 }
292 }
293
294
point_binning_allocate(struct PointBinning * point_binning,int rows,int cols,RASTER_MAP_TYPE rtype)295 void point_binning_allocate(struct PointBinning *point_binning, int rows,
296 int cols, RASTER_MAP_TYPE rtype)
297 {
298 if (point_binning->bin_n) {
299 G_debug(2, "allocating n_array");
300 point_binning->n_array =
301 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(CELL_TYPE));
302 blank_array(point_binning->n_array, rows, cols, CELL_TYPE, 0);
303 }
304 if (point_binning->bin_min) {
305 G_debug(2, "allocating min_array");
306 point_binning->min_array =
307 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
308 blank_array(point_binning->min_array, rows, cols, rtype, -1); /* fill with NULLs */
309 }
310 if (point_binning->bin_max) {
311 G_debug(2, "allocating max_array");
312 point_binning->max_array =
313 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
314 blank_array(point_binning->max_array, rows, cols, rtype, -1); /* fill with NULLs */
315 }
316 if (point_binning->bin_sum) {
317 G_debug(2, "allocating sum_array");
318 point_binning->sum_array =
319 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
320 blank_array(point_binning->sum_array, rows, cols, rtype, 0);
321 }
322 if (point_binning->bin_sumsq) {
323 G_debug(2, "allocating sumsq_array");
324 point_binning->sumsq_array =
325 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
326 blank_array(point_binning->sumsq_array, rows, cols, rtype, 0);
327 }
328 if (point_binning->bin_index) {
329 G_debug(2, "allocating index_array");
330 point_binning->index_array =
331 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(CELL_TYPE));
332 blank_array(point_binning->index_array, rows, cols, CELL_TYPE, -1); /* fill with NULLs */
333 }
334 if (point_binning->bin_coordinates) {
335 G_debug(2, "allocating x_array and y_array");
336 point_binning->x_array =
337 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
338 blank_array(point_binning->x_array, rows, cols, rtype, 0);
339 point_binning->y_array =
340 G_calloc((size_t) rows * (cols + 1), Rast_cell_size(rtype));
341 blank_array(point_binning->y_array, rows, cols, rtype, 0);
342 }
343 }
344
point_binning_free(struct PointBinning * point_binning,struct BinIndex * bin_index_nodes)345 void point_binning_free(struct PointBinning *point_binning,
346 struct BinIndex *bin_index_nodes)
347 {
348 if (point_binning->bin_n)
349 G_free(point_binning->n_array);
350 if (point_binning->bin_min)
351 G_free(point_binning->min_array);
352 if (point_binning->bin_max)
353 G_free(point_binning->max_array);
354 if (point_binning->bin_sum)
355 G_free(point_binning->sum_array);
356 if (point_binning->bin_sumsq)
357 G_free(point_binning->sumsq_array);
358 if (point_binning->bin_index) {
359 G_free(point_binning->index_array);
360 G_free(bin_index_nodes->nodes);
361 bin_index_nodes->num_nodes = 0;
362 bin_index_nodes->max_nodes = 0;
363 bin_index_nodes->nodes = NULL;
364 }
365 if (point_binning->bin_coordinates) {
366 G_free(point_binning->x_array);
367 G_free(point_binning->y_array);
368 }
369 }
370
write_variance(void * raster_row,void * n_array,void * sum_array,void * sumsq_array,int row,int cols,RASTER_MAP_TYPE rtype,int method)371 void write_variance(void *raster_row, void *n_array, void *sum_array,
372 void *sumsq_array, int row, int cols,
373 RASTER_MAP_TYPE rtype, int method)
374 {
375 size_t offset, n_offset;
376 int n = 0;
377 double variance;
378 double sum = 0.;
379 double sumsq = 0.;
380 int col;
381 void *ptr = raster_row;
382
383 for (col = 0; col < cols; col++) {
384 offset = ((size_t) row * cols + col) * Rast_cell_size(rtype);
385 n_offset = ((size_t) row * cols + col) * Rast_cell_size(CELL_TYPE);
386 n = Rast_get_c_value(n_array + n_offset, CELL_TYPE);
387 sum = Rast_get_d_value(sum_array + offset, rtype);
388 sumsq = Rast_get_d_value(sumsq_array + offset, rtype);
389
390 if (n == 0)
391 Rast_set_null_value(ptr, 1, rtype);
392 else if (n == 1)
393 Rast_set_d_value(ptr, 0.0, rtype);
394 else {
395 variance = (sumsq - sum * sum / n) / n;
396 if (variance < GRASS_EPSILON)
397 variance = 0.0;
398
399 /* nan test */
400 if (variance != variance)
401 Rast_set_null_value(ptr, 1, rtype);
402 else {
403
404 if (method == METHOD_STDDEV)
405 variance = sqrt(variance);
406
407 else if (method == METHOD_COEFF_VAR)
408 variance = 100 * sqrt(variance) / (sum / n);
409
410 /* nan test */
411 if (variance != variance)
412 variance = 0.0; /* OK for n > 0 ? */
413
414 Rast_set_d_value(ptr, variance, rtype);
415 }
416
417 }
418 ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
419 }
420 }
421
write_median(struct BinIndex * bin_index,void * raster_row,void * index_array,int row,int cols,RASTER_MAP_TYPE rtype)422 void write_median(struct BinIndex *bin_index, void *raster_row,
423 void *index_array, int row, int cols, RASTER_MAP_TYPE rtype)
424 {
425 size_t n_offset;
426 int n;
427 int j;
428 double z;
429 int col;
430 int node_id, head_id;
431 void *ptr = raster_row;
432
433 for (col = 0; col < cols; col++) {
434 n_offset = ((size_t) row * cols + col) * Rast_cell_size(CELL_TYPE);
435 if (Rast_is_null_value(index_array + n_offset, CELL_TYPE)) /* no points in cell */
436 Rast_set_null_value(ptr, 1, rtype);
437 else { /* one or more points in cell */
438
439 head_id = Rast_get_c_value(index_array + n_offset, CELL_TYPE);
440 node_id = head_id;
441
442 n = 0;
443
444 while (node_id != -1) { /* count number of points in cell */
445 n++;
446 node_id = bin_index->nodes[node_id].next;
447 }
448
449 if (n == 1) /* only one point, use that */
450 Rast_set_d_value(ptr, bin_index->nodes[head_id].z, rtype);
451 else if (n % 2 != 0) { /* odd number of points: median_i = (n + 1) / 2 */
452 n = (n + 1) / 2;
453 node_id = head_id;
454 for (j = 1; j < n; j++) /* get "median element" */
455 node_id = bin_index->nodes[node_id].next;
456
457 Rast_set_d_value(ptr, bin_index->nodes[node_id].z, rtype);
458 }
459 else { /* even number of points: median = (val_below + val_above) / 2 */
460
461 z = (n + 1) / 2.0;
462 n = floor(z);
463 node_id = head_id;
464 for (j = 1; j < n; j++) /* get element "below" */
465 node_id = bin_index->nodes[node_id].next;
466
467 z = (bin_index->nodes[node_id].z +
468 bin_index->nodes[bin_index->nodes[node_id].next].z) / 2;
469 Rast_set_d_value(ptr, z, rtype);
470 }
471 }
472 ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
473 }
474 }
475
write_percentile(struct BinIndex * bin_index,void * raster_row,void * index_array,int row,int cols,RASTER_MAP_TYPE rtype,int pth)476 void write_percentile(struct BinIndex *bin_index, void *raster_row,
477 void *index_array, int row, int cols,
478 RASTER_MAP_TYPE rtype, int pth)
479 {
480 size_t n_offset;
481 int n;
482 int j;
483 double z;
484 int col;
485 int node_id, head_id;
486 int r_low, r_up;
487 void *ptr = raster_row;
488
489 for (col = 0; col < cols; col++) {
490 n_offset = ((size_t) row * cols + col) * Rast_cell_size(CELL_TYPE);
491 if (Rast_is_null_value(index_array + n_offset, CELL_TYPE)) /* no points in cell */
492 Rast_set_null_value(ptr, 1, rtype);
493 else {
494 head_id = Rast_get_c_value(index_array + n_offset, CELL_TYPE);
495 node_id = head_id;
496 n = 0;
497
498 while (node_id != -1) { /* count number of points in cell */
499 n++;
500 node_id = bin_index->nodes[node_id].next;
501 }
502
503 z = (pth * (n + 1)) / 100.0;
504 r_low = floor(z); /* lower rank */
505 if (r_low < 1)
506 r_low = 1;
507 else if (r_low > n)
508 r_low = n;
509
510 r_up = ceil(z); /* upper rank */
511 if (r_up > n)
512 r_up = n;
513
514 node_id = head_id;
515 for (j = 1; j < r_low; j++) /* search lower value */
516 node_id = bin_index->nodes[node_id].next;
517
518 z = bin_index->nodes[node_id].z; /* save lower value */
519 node_id = head_id;
520 for (j = 1; j < r_up; j++) /* search upper value */
521 node_id = bin_index->nodes[node_id].next;
522
523 z = (z + bin_index->nodes[node_id].z) / 2;
524 Rast_set_d_value(ptr, z, rtype);
525 }
526 ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
527 }
528 }
529
write_skewness(struct BinIndex * bin_index,void * raster_row,void * index_array,int row,int cols,RASTER_MAP_TYPE rtype)530 void write_skewness(struct BinIndex *bin_index, void *raster_row,
531 void *index_array, int row, int cols,
532 RASTER_MAP_TYPE rtype)
533 {
534 size_t n_offset;
535 int n;
536 double z;
537 int col;
538 int node_id, head_id;
539 double variance, mean, skew, sumdev;
540 double sum = 0.;
541 double sumsq = 0.;
542 void *ptr = raster_row;
543
544 for (col = 0; col < cols; col++) {
545 n_offset = ((size_t) row * cols + col) * Rast_cell_size(CELL_TYPE);
546 if (Rast_is_null_value(index_array + n_offset, CELL_TYPE)) /* no points in cell */
547 Rast_set_null_value(ptr, 1, rtype);
548 else {
549 head_id = Rast_get_c_value(index_array + n_offset, CELL_TYPE);
550 node_id = head_id;
551
552 n = 0; /* count */
553 sum = 0.0; /* sum */
554 sumsq = 0.0; /* sum of squares */
555 sumdev = 0.0; /* sum of (xi - mean)^3 */
556 skew = 0.0; /* skewness */
557
558 while (node_id != -1) {
559 z = bin_index->nodes[node_id].z;
560 n++;
561 sum += z;
562 sumsq += (z * z);
563 node_id = bin_index->nodes[node_id].next;
564 }
565
566 if (n > 1) { /* if n == 1, skew is "0.0" */
567 mean = sum / n;
568 node_id = head_id;
569 while (node_id != -1) {
570 z = bin_index->nodes[node_id].z;
571 sumdev += pow((z - mean), 3);
572 node_id = bin_index->nodes[node_id].next;
573 }
574
575 variance = (sumsq - sum * sum / n) / n;
576 if (variance < GRASS_EPSILON)
577 skew = 0.0;
578 else
579 skew = sumdev / ((n - 1) * pow(sqrt(variance), 3));
580 }
581 Rast_set_d_value(ptr, skew, rtype);
582 }
583 ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
584 }
585 }
586
write_trimmean(struct BinIndex * bin_index,void * raster_row,void * index_array,int row,int cols,RASTER_MAP_TYPE rtype,double trim)587 void write_trimmean(struct BinIndex *bin_index, void *raster_row,
588 void *index_array, int row, int cols,
589 RASTER_MAP_TYPE rtype, double trim)
590 {
591 size_t n_offset;
592 int n;
593 int j, k;
594 int col;
595 int node_id, head_id;
596 double mean;
597 double sum = 0.;
598 void *ptr = raster_row;
599
600 for (col = 0; col < cols; col++) {
601 n_offset = ((size_t) row * cols + col) * Rast_cell_size(CELL_TYPE);
602 if (Rast_is_null_value(index_array + n_offset, CELL_TYPE)) /* no points in cell */
603 Rast_set_null_value(ptr, 1, rtype);
604 else {
605 head_id = Rast_get_c_value(index_array + n_offset, CELL_TYPE);
606
607 node_id = head_id;
608 n = 0;
609 while (node_id != -1) { /* count number of points in cell */
610 n++;
611 node_id = bin_index->nodes[node_id].next;
612 }
613
614 if (1 == n)
615 mean = bin_index->nodes[head_id].z;
616 else {
617 k = floor(trim * n + 0.5); /* number of ranks to discard on each tail */
618
619 if (k > 0 && (n - 2 * k) > 0) { /* enough elements to discard */
620 node_id = head_id;
621 for (j = 0; j < k; j++) /* move to first rank to consider */
622 node_id = bin_index->nodes[node_id].next;
623
624 j = k + 1;
625 k = n - k;
626 n = 0;
627 sum = 0.0;
628
629 while (j <= k) { /* get values in interval */
630 n++;
631 sum += bin_index->nodes[node_id].z;
632 node_id = bin_index->nodes[node_id].next;
633 j++;
634 }
635 }
636 else {
637 node_id = head_id;
638 n = 0;
639 sum = 0.0;
640 while (node_id != -1) {
641 n++;
642 sum += bin_index->nodes[node_id].z;
643 node_id = bin_index->nodes[node_id].next;
644 }
645 }
646 mean = sum / n;
647 }
648 Rast_set_d_value(ptr, mean, rtype);
649 }
650 ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
651 }
652 }
653
write_values(struct PointBinning * point_binning,struct BinIndex * bin_index_nodes,void * raster_row,int row,int cols,RASTER_MAP_TYPE rtype,struct VectorWriter * vector_writer)654 void write_values(struct PointBinning *point_binning,
655 struct BinIndex *bin_index_nodes, void *raster_row, int row,
656 int cols, RASTER_MAP_TYPE rtype,
657 struct VectorWriter *vector_writer)
658 {
659 void *ptr = NULL;
660 int col;
661
662 switch (point_binning->method) {
663 case METHOD_N: /* n is a straight copy */
664 Rast_raster_cpy(raster_row,
665 point_binning->n_array +
666 ((size_t) row * cols * Rast_cell_size(CELL_TYPE)), cols,
667 CELL_TYPE);
668 break;
669
670 case METHOD_MIN:
671 Rast_raster_cpy(raster_row,
672 point_binning->min_array +
673 ((size_t) row * cols * Rast_cell_size(rtype)), cols, rtype);
674 break;
675
676 case METHOD_MAX:
677 Rast_raster_cpy(raster_row,
678 point_binning->max_array +
679 ((size_t) row * cols * Rast_cell_size(rtype)), cols, rtype);
680 break;
681
682 case METHOD_SUM:
683 Rast_raster_cpy(raster_row,
684 point_binning->sum_array +
685 ((size_t) row * cols * Rast_cell_size(rtype)), cols, rtype);
686 break;
687
688 case METHOD_RANGE: /* (max-min) */
689 ptr = raster_row;
690 for (col = 0; col < cols; col++) {
691 size_t offset = ((size_t) row * cols + col) * Rast_cell_size(rtype);
692 double min =
693 Rast_get_d_value(point_binning->min_array + offset, rtype);
694 double max =
695 Rast_get_d_value(point_binning->max_array + offset, rtype);
696 Rast_set_d_value(ptr, max - min, rtype);
697 ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
698 }
699 break;
700
701 case METHOD_MEAN: /* (sum / n) */
702 ptr = raster_row;
703 for (col = 0; col < cols; col++) {
704 size_t offset = ((size_t) row * cols + col) * Rast_cell_size(rtype);
705 size_t n_offset = ((size_t) row * cols + col) * Rast_cell_size(CELL_TYPE);
706 int n = Rast_get_c_value(point_binning->n_array + n_offset,
707 CELL_TYPE);
708 double sum =
709 Rast_get_d_value(point_binning->sum_array + offset, rtype);
710
711 if (n == 0)
712 Rast_set_null_value(ptr, 1, rtype);
713 else
714 Rast_set_d_value(ptr, (sum / n), rtype);
715
716 ptr = G_incr_void_ptr(ptr, Rast_cell_size(rtype));
717 }
718 break;
719
720 case METHOD_STDDEV: /* sqrt(variance) */
721 case METHOD_VARIANCE: /* (sumsq - sum*sum/n)/n */
722 case METHOD_COEFF_VAR: /* 100 * stdev / mean */
723 write_variance(raster_row, point_binning->n_array,
724 point_binning->sum_array, point_binning->sumsq_array,
725 row, cols, rtype, point_binning->method);
726 break;
727 case METHOD_MEDIAN: /* median, if only one point in cell we will use that */
728 write_median(bin_index_nodes, raster_row, point_binning->index_array,
729 row, cols, rtype);
730 break;
731 case METHOD_PERCENTILE: /* rank = (pth*(n+1))/100; interpolate linearly */
732 write_percentile(bin_index_nodes, raster_row,
733 point_binning->index_array, row, cols, rtype,
734 point_binning->pth);
735 break;
736 case METHOD_SKEWNESS: /* skewness = sum(xi-mean)^3/(N-1)*s^3 */
737 write_skewness(bin_index_nodes, raster_row,
738 point_binning->index_array, row, cols, rtype);
739 break;
740 case METHOD_TRIMMEAN:
741 write_trimmean(bin_index_nodes, raster_row,
742 point_binning->index_array, row, cols, rtype,
743 point_binning->trim);
744 break;
745
746 default:
747 G_debug(2, "No method selected");
748 }
749 if (point_binning->bin_coordinates) {
750 for (col = 0; col < cols; col++) {
751 size_t offset = ((size_t) row * cols + col) * Rast_cell_size(rtype);
752 size_t n_offset = ((size_t) row * cols + col) * Rast_cell_size(CELL_TYPE);
753 int n = Rast_get_c_value(point_binning->n_array + n_offset,
754 CELL_TYPE);
755
756 if (n == 0)
757 continue;
758
759 double sum_x =
760 Rast_get_d_value(point_binning->x_array + offset, rtype);
761 double sum_y =
762 Rast_get_d_value(point_binning->y_array + offset, rtype);
763 /* TODO: we do this also in mean writing */
764 double sum_z =
765 Rast_get_d_value(point_binning->sum_array + offset, rtype);
766
767 /* We are not writing any categories. They are not needed
768 * and potentially it is too much trouble to do it and it is
769 * unclear what to write. Not writing them is also a little
770 * bit faster. */
771 Vect_append_point(vector_writer->points, sum_x, sum_y, sum_z / n);
772 Vect_write_line(vector_writer->info, GV_POINT,
773 vector_writer->points, vector_writer->cats);
774 Vect_reset_line(vector_writer->points);
775 vector_writer->count++;
776 }
777 }
778 }
779
780 /* TODO: duplication with support.c, refactoring needed */
get_cell_ptr(void * array,int cols,int row,int col,RASTER_MAP_TYPE map_type)781 static void *get_cell_ptr(void *array, int cols, int row, int col,
782 RASTER_MAP_TYPE map_type)
783 {
784 return G_incr_void_ptr(array,
785 ((row * (size_t) cols) +
786 col) * Rast_cell_size(map_type));
787 }
788
update_val(void * array,int cols,int row,int col,RASTER_MAP_TYPE map_type,double value)789 int update_val(void *array, int cols, int row, int col,
790 RASTER_MAP_TYPE map_type, double value)
791 {
792 void *ptr = get_cell_ptr(array, cols, row, col, map_type);
793
794 Rast_set_d_value(ptr, value, map_type);
795 return 0;
796 }
797
update_moving_mean(void * array,int cols,int row,int col,RASTER_MAP_TYPE rtype,double value,int n)798 int update_moving_mean(void *array, int cols, int row, int col,
799 RASTER_MAP_TYPE rtype, double value, int n)
800 {
801 /* for xy we do this check twice */
802 if (n != 0) {
803 double m_v;
804
805 row_array_get_value_row_col(array, row, col, cols, rtype, &m_v);
806 value = m_v + (value - m_v) / n;
807 }
808 /* else we just write the initial value */
809 return update_val(array, cols, row, col, rtype, value);;
810 }
811
update_value(struct PointBinning * point_binning,struct BinIndex * bin_index_nodes,int cols,int arr_row,int arr_col,RASTER_MAP_TYPE rtype,double x,double y,double z)812 void update_value(struct PointBinning *point_binning,
813 struct BinIndex *bin_index_nodes, int cols, int arr_row,
814 int arr_col, RASTER_MAP_TYPE rtype, double x, double y,
815 double z)
816 {
817 if (point_binning->bin_n)
818 update_n(point_binning->n_array, cols, arr_row, arr_col);
819 if (point_binning->bin_min)
820 update_min(point_binning->min_array, cols, arr_row, arr_col, rtype,
821 z);
822 if (point_binning->bin_max)
823 update_max(point_binning->max_array, cols, arr_row, arr_col, rtype,
824 z);
825 if (point_binning->bin_sum)
826 update_sum(point_binning->sum_array, cols, arr_row, arr_col, rtype,
827 z);
828 if (point_binning->bin_sumsq)
829 update_sumsq(point_binning->sumsq_array, cols, arr_row, arr_col,
830 rtype, z);
831 if (point_binning->bin_index)
832 update_bin_index(bin_index_nodes, point_binning->index_array, cols,
833 arr_row, arr_col, rtype, z);
834 if (point_binning->bin_coordinates) {
835 /* this assumes that n is already computed for this xyz */
836 void *ptr = get_cell_ptr(point_binning->n_array, cols, arr_row,
837 arr_col, CELL_TYPE);
838 int n = Rast_get_c_value(ptr, CELL_TYPE);
839
840 update_moving_mean(point_binning->x_array, cols, arr_row, arr_col,
841 rtype, x, n);
842 update_moving_mean(point_binning->y_array, cols, arr_row, arr_col,
843 rtype, y, n);
844 }
845 }
846