GRASS GIS 8 Programmer's Manual 8.2.0(2022)-exported
n_arrays_calc.c
Go to the documentation of this file.
1
2/*****************************************************************************
3*
4* MODULE: Grass PDE Numerical Library
5* AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6* soerengebbert <at> gmx <dot> de
7*
8* PURPOSE: Higher level array management functions
9* part of the gpde library
10*
11* COPYRIGHT: (C) 2000 by the GRASS Development Team
12*
13* This program is free software under the GNU General Public
14* License (>=v2). Read the file COPYING that comes with GRASS
15* for details.
16*
17*****************************************************************************/
18#include <math.h>
19
20#include <grass/N_pde.h>
21#include <grass/raster.h>
22#include <grass/glocale.h>
23
24
25/* ******************** 2D ARRAY FUNCTIONS *********************** */
26
27/*!
28 * \brief Copy the source N_array_2d struct to the target N_array_2d struct
29 *
30 * The arrays must have the same size and the same offset.
31 *
32 * The array types can be mixed, the values are automatically casted
33 * and the null values are set accordingly.
34 * <br><br>
35 * If you copy a cell array into a dcell array, the values are casted to dcell and
36 * the null values are converted from cell-null to dcell-null
37 * <br><br>
38 * This function can be called in a parallel region defined with OpenMP.
39 * The copy loop is parallelize with a openmp for pragma.
40 *
41 * \param source N_array_2d *
42 * \param target N_array_2d *
43 * \return void
44 * */
46{
47 int i;
48 int null = 0;
49
50#pragma omp single
51 {
52 if (source->cols_intern != target->cols_intern)
54 ("N_copy_array_2d: the arrays are not of equal size");
55
56 if (source->rows_intern != target->rows_intern)
58 ("N_copy_array_2d: the arrays are not of equal size");
59
60 G_debug(3,
61 "N_copy_array_2d: copy source array to target array size %i",
62 source->cols_intern * source->rows_intern);
63 }
64
65#pragma omp for
66 for (i = 0; i < source->cols_intern * source->rows_intern; i++) {
67 null = 0;
68 if (source->type == CELL_TYPE) {
69 if (Rast_is_c_null_value((void *)&source->cell_array[i]))
70 null = 1;
71
72 if (target->type == CELL_TYPE) {
73 target->cell_array[i] = source->cell_array[i];
74 }
75 if (target->type == FCELL_TYPE) {
76 if (null)
77 Rast_set_f_null_value((void *)&(target->fcell_array[i]), 1);
78 else
79 target->fcell_array[i] = (FCELL) source->cell_array[i];
80 }
81 if (target->type == DCELL_TYPE) {
82 if (null)
83 Rast_set_d_null_value((void *)&(target->dcell_array[i]), 1);
84 else
85 target->dcell_array[i] = (DCELL) source->cell_array[i];
86 }
87
88 }
89 if (source->type == FCELL_TYPE) {
90 if (Rast_is_f_null_value((void *)&source->fcell_array[i]))
91 null = 1;
92
93 if (target->type == CELL_TYPE) {
94 if (null)
95 Rast_set_c_null_value((void *)&(target->cell_array[i]), 1);
96 else
97 target->cell_array[i] = (CELL) source->fcell_array[i];
98 }
99 if (target->type == FCELL_TYPE) {
100 target->fcell_array[i] = source->fcell_array[i];
101 }
102 if (target->type == DCELL_TYPE) {
103 if (null)
104 Rast_set_d_null_value((void *)&(target->dcell_array[i]), 1);
105 else
106 target->dcell_array[i] = (DCELL) source->fcell_array[i];
107 }
108 }
109 if (source->type == DCELL_TYPE) {
110 if (Rast_is_d_null_value((void *)&source->dcell_array[i]))
111 null = 1;
112
113 if (target->type == CELL_TYPE) {
114 if (null)
115 Rast_set_c_null_value((void *)&(target->cell_array[i]), 1);
116 else
117 target->cell_array[i] = (CELL) source->dcell_array[i];
118 }
119 if (target->type == FCELL_TYPE) {
120 if (null)
121 Rast_set_f_null_value((void *)&(target->fcell_array[i]), 1);
122 else
123 target->fcell_array[i] = (FCELL) source->dcell_array[i];
124 }
125 if (target->type == DCELL_TYPE) {
126 target->dcell_array[i] = source->dcell_array[i];
127 }
128 }
129 }
130
131 return;
132}
133
134/*!
135 * \brief Calculate the norm of the two input arrays
136 *
137 * The norm can be of type N_MAXIMUM_NORM or N_EUKLID_NORM.
138 * All arrays must have equal sizes and offsets.
139 * The complete data array inclusively offsets is used for norm calucaltion.
140 * Only non-null values are used to calculate the norm.
141 *
142
143 * \param a N_array_2d *
144 * \param b N_array_2d *
145 * \param type the type of the norm -> N_MAXIMUM_NORM, N_EUKLID_NORM
146 * \return double the calculated norm
147 * */
148double N_norm_array_2d(N_array_2d * a, N_array_2d * b, int type)
149{
150 int i = 0;
151 double norm = 0.0, tmp = 0.0;
152 double v1 = 0.0, v2 = 0.0;
153
154 if (a->cols_intern != b->cols_intern)
155 G_fatal_error("N_norm_array_2d: the arrays are not of equal size");
156
157 if (a->rows_intern != b->rows_intern)
158 G_fatal_error("N_norm_array_2d: the arrays are not of equal size");
159
160 G_debug(3, "N_norm_array_2d: norm of a and b size %i",
161 a->cols_intern * a->rows_intern);
162
163 for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
164 v1 = 0.0;
165 v2 = 0.0;
166
167 if (a->type == CELL_TYPE) {
168 if (!Rast_is_f_null_value((void *)&(a->cell_array[i])))
169 v1 = (double)a->cell_array[i];
170 }
171 if (a->type == FCELL_TYPE) {
172 if (!Rast_is_f_null_value((void *)&(a->fcell_array[i])))
173 v1 = (double)a->fcell_array[i];
174 }
175 if (a->type == DCELL_TYPE) {
176 if (!Rast_is_f_null_value((void *)&(a->dcell_array[i])))
177 v1 = (double)a->dcell_array[i];
178 }
179 if (b->type == CELL_TYPE) {
180 if (!Rast_is_f_null_value((void *)&(b->cell_array[i])))
181 v2 = (double)b->cell_array[i];
182 }
183 if (b->type == FCELL_TYPE) {
184 if (!Rast_is_f_null_value((void *)&(b->fcell_array[i])))
185 v2 = (double)b->fcell_array[i];
186 }
187 if (b->type == DCELL_TYPE) {
188 if (!Rast_is_f_null_value((void *)&(b->dcell_array[i])))
189 v2 = (double)b->dcell_array[i];
190 }
191
192 if (type == N_MAXIMUM_NORM) {
193 tmp = fabs(v2 - v1);
194 if ((tmp > norm))
195 norm = tmp;
196 }
197 if (type == N_EUKLID_NORM) {
198 norm += fabs(v2 - v1);
199 }
200 }
201
202 return norm;
203}
204
205/*!
206 * \brief Calculate basic statistics of the N_array_2d struct
207 *
208 * Calculates the minimum, maximum, sum and the number of
209 * non null values. The array offset can be included in the calculation.
210 *
211 * \param a N_array_2d * - input array
212 * \param min double* - variable to store the computed minimum
213 * \param max double* - variable to store the computed maximum
214 * \param sum double* - variable to store the computed sum
215 * \param nonull int* - variable to store the number of non null values
216 * \param withoffset - if 1 include offset values in statistic calculation, 0 otherwise
217 * \return void
218 * */
219void N_calc_array_2d_stats(N_array_2d * a, double *min, double *max,
220 double *sum, int *nonull, int withoffset)
221{
222 int i, j;
223 double val;
224
225 *sum = 0.0;
226 *nonull = 0;
227
228 if (withoffset == 1) {
229
230 *min =
231 (double)N_get_array_2d_d_value(a, 0 - a->offset, 0 - a->offset);
232 *max =
233 (double)N_get_array_2d_d_value(a, 0 - a->offset, 0 - a->offset);
234
235 for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
236 for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
237 if (!N_is_array_2d_value_null(a, i, j)) {
238 val = (double)N_get_array_2d_d_value(a, i, j);
239 if (*min > val)
240 *min = val;
241 if (*max < val)
242 *max = val;
243 *sum += val;
244 (*nonull)++;
245 }
246 }
247 }
248 }
249 else {
250
251 *min = (double)N_get_array_2d_d_value(a, 0, 0);
252 *max = (double)N_get_array_2d_d_value(a, 0, 0);
253
254
255 for (j = 0; j < a->rows; j++) {
256 for (i = 0; i < a->cols; i++) {
257 if (!N_is_array_2d_value_null(a, i, j)) {
258 val = (double)N_get_array_2d_d_value(a, i, j);
259 if (*min > val)
260 *min = val;
261 if (*max < val)
262 *max = val;
263 *sum += val;
264 (*nonull)++;
265 }
266 }
267 }
268 }
269
270 G_debug(3,
271 "N_calc_array_2d_stats: compute array stats, min %g, max %g, sum %g, nonull %i",
272 *min, *max, *sum, *nonull);
273 return;
274}
275
276
277/*!
278 * \brief Perform calculations with two input arrays,
279 * the result is written to a third array.
280 *
281 * All arrays must have equal sizes and offsets.
282 * The complete data array inclusively offsets is used for calucaltions.
283 * Only non-null values are computed. If one array value is null,
284 * the result array value will be null too.
285 * <br><br>
286 * If a division with zero is detected, the resulting arrays
287 * value will set to null and not to NaN.
288 * <br><br>
289 * The result array is optional, if the result arrays points to NULL,
290 * a new array will be allocated with the largest arrays data type
291 * (CELL, FCELL or DCELL) used by the input arrays.
292 * <br><br>
293 * the array computations can be of the following forms:
294 *
295 * <ul>
296 * <li>result = a + b -> N_ARRAY_SUM</li>
297 * <li>result = a - b -> N_ARRAY_DIF</li>
298 * <li>result = a * b -> N_ARRAY_MUL</li>
299 * <li>result = a / b -> N_ARRAY_DIV</li>
300 * </ul>
301 *
302 * \param a N_array_2d * - first input array
303 * \param b N_array_2d * - second input array
304 * \param result N_array_2d * - the optional result array
305 * \param type - the type of calculation
306 * \return N_array_2d * - the pointer to the result array
307 * */
309 N_array_2d * result, int type)
310{
311 N_array_2d *c;
312 int i, j, setnull = 0;
313 double va = 0.0, vb = 0.0, vc = 0.0; /*variables used for calculation */
314
315 /*Set the pointer */
316 c = result;
317
318#pragma omp single
319 {
320 /*Check the array sizes */
321 if (a->cols_intern != b->cols_intern)
323 ("N_math_array_2d: the arrays are not of equal size");
324 if (a->rows_intern != b->rows_intern)
326 ("N_math_array_2d: the arrays are not of equal size");
327 if (a->offset != b->offset)
329 ("N_math_array_2d: the arrays have different offsets");
330
331 G_debug(3, "N_math_array_2d: mathematical calculations, size: %i",
332 a->cols_intern * a->rows_intern);
333
334 /*if the result array is null, allocate a new one, use the
335 * largest data type of the input arrays*/
336 if (c == NULL) {
337 if (a->type == DCELL_TYPE || b->type == DCELL_TYPE) {
338 c = N_alloc_array_2d(a->cols, a->rows, a->offset, DCELL_TYPE);
339 G_debug(3,
340 "N_math_array_2d: array of type DCELL_TYPE created");
341 }
342 else if (a->type == FCELL_TYPE || b->type == FCELL_TYPE) {
343 c = N_alloc_array_2d(a->cols, a->rows, a->offset, FCELL_TYPE);
344 G_debug(3,
345 "N_math_array_2d: array of type FCELL_TYPE created");
346 }
347 else {
348 c = N_alloc_array_2d(a->cols, a->rows, a->offset, CELL_TYPE);
349 G_debug(3,
350 "N_math_array_2d: array of type CELL_TYPE created");
351 }
352 }
353 else {
354 /*Check the array sizes */
355 if (a->cols_intern != c->cols_intern)
357 ("N_math_array_2d: the arrays are not of equal size");
358 if (a->rows_intern != c->rows_intern)
360 ("N_math_array_2d: the arrays are not of equal size");
361 if (a->offset != c->offset)
363 ("N_math_array_2d: the arrays have different offsets");
364 }
365 }
366
367#pragma omp for private(va, vb, vc, setnull)
368 for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
369 for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
370 if (!N_is_array_2d_value_null(a, i, j) &&
371 !N_is_array_2d_value_null(b, i, j)) {
372 /*we always calculate internally with double values */
373 va = (double)N_get_array_2d_d_value(a, i, j);
374 vb = (double)N_get_array_2d_d_value(b, i, j);
375 vc = 0;
376 setnull = 0;
377
378 switch (type) {
379 case N_ARRAY_SUM:
380 vc = va + vb;
381 break;
382 case N_ARRAY_DIF:
383 vc = va - vb;
384 break;
385 case N_ARRAY_MUL:
386 vc = va * vb;
387 break;
388 case N_ARRAY_DIV:
389 if (vb != 0)
390 vc = va / vb;
391 else
392 setnull = 1;
393 break;
394 }
395
396 if (c->type == CELL_TYPE) {
397 if (setnull)
399 else
400 N_put_array_2d_c_value(c, i, j, (CELL) vc);
401 }
402 if (c->type == FCELL_TYPE) {
403 if (setnull)
405 else
406 N_put_array_2d_f_value(c, i, j, (FCELL) vc);
407 }
408 if (c->type == DCELL_TYPE) {
409 if (setnull)
411 else
412 N_put_array_2d_d_value(c, i, j, (DCELL) vc);
413 }
414
415 }
416 else {
418 }
419 }
420 }
421
422 return c;
423}
424
425/*!
426 * \brief Convert all null values to zero values
427 *
428 * The complete data array inclusively offsets is used.
429 * The array data types are automatically recognized.
430 *
431 * \param a N_array_2d *
432 * \return int - number of replaced values
433 * */
435{
436 int i = 0, count = 0;
437
438 G_debug(3, "N_convert_array_2d_null_to_zero: convert array of size %i",
439 a->cols_intern * a->rows_intern);
440
441 if (a->type == CELL_TYPE)
442 for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
443 if (Rast_is_c_null_value((void *)&(a->cell_array[i]))) {
444 a->cell_array[i] = 0;
445 count++;
446 }
447 }
448
449 if (a->type == FCELL_TYPE)
450 for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
451 if (Rast_is_f_null_value((void *)&(a->fcell_array[i]))) {
452 a->fcell_array[i] = 0.0;
453 count++;
454 }
455 }
456
457
458 if (a->type == DCELL_TYPE)
459 for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
460 if (Rast_is_d_null_value((void *)&(a->dcell_array[i]))) {
461 a->dcell_array[i] = 0.0;
462 count++;
463 }
464 }
465
466
467 if (a->type == CELL_TYPE)
468 G_debug(2,
469 "N_convert_array_2d_null_to_zero: %i values of type CELL_TYPE are converted",
470 count);
471 if (a->type == FCELL_TYPE)
472 G_debug(2,
473 "N_convert_array_2d_null_to_zero: %i valuess of type FCELL_TYPE are converted",
474 count);
475 if (a->type == DCELL_TYPE)
476 G_debug(2,
477 "N_convert_array_2d_null_to_zero: %i valuess of type DCELL_TYPE are converted",
478 count);
479
480 return count;
481}
482
483/* ******************** 3D ARRAY FUNCTIONS *********************** */
484
485/*!
486 * \brief Copy the source N_array_3d struct to the target N_array_3d struct
487 *
488 * The arrays must have the same size and the same offset.
489 *
490 * The array data types can be mixed, the values are automatically casted
491 * and the null values are set accordingly.
492 *
493 * If you copy a float array to a double array, the values are casted to DCELL and
494 * the null values are converted from FCELL-null to DCELL-null
495 *
496 * \param source N_array_3d *
497 * \param target N_array_3d *
498 * \return void
499 * */
501{
502 int i;
503 int null;
504
505 if (source->cols_intern != target->cols_intern)
506 G_fatal_error("N_copy_array_3d: the arrays are not of equal size");
507
508 if (source->rows_intern != target->rows_intern)
509 G_fatal_error("N_copy_array_3d: the arrays are not of equal size");
510
511 if (source->depths_intern != target->depths_intern)
512 G_fatal_error("N_copy_array_3d: the arrays are not of equal size");
513
514
515 G_debug(3, "N_copy_array_3d: copy source array to target array size %i",
516 source->cols_intern * source->rows_intern *
517 source->depths_intern);
518
519 for (i = 0;
520 i <
521 source->cols_intern * source->rows_intern * source->depths_intern;
522 i++) {
523 null = 0;
524 if (source->type == FCELL_TYPE) {
525 if (Rast3d_is_null_value_num
526 ((void *)&(source->fcell_array[i]), FCELL_TYPE))
527 null = 1;
528
529 if (target->type == FCELL_TYPE) {
530 target->fcell_array[i] = source->fcell_array[i];
531 }
532 if (target->type == DCELL_TYPE) {
533 if (null)
534 Rast3d_set_null_value((void *)&(target->dcell_array[i]), 1,
535 DCELL_TYPE);
536 else
537 target->dcell_array[i] = (double)source->fcell_array[i];
538 }
539
540 }
541 if (source->type == DCELL_TYPE) {
542 if (Rast3d_is_null_value_num
543 ((void *)&(source->dcell_array[i]), DCELL_TYPE))
544 null = 1;
545
546 if (target->type == FCELL_TYPE) {
547 if (null)
548 Rast3d_set_null_value((void *)&(target->fcell_array[i]), 1,
549 FCELL_TYPE);
550 else
551 target->fcell_array[i] = (float)source->dcell_array[i];
552 }
553 if (target->type == DCELL_TYPE) {
554 target->dcell_array[i] = source->dcell_array[i];
555 }
556 }
557 }
558
559 return;
560}
561
562
563/*!
564 * \brief Calculate the norm of the two input arrays
565 *
566 * The norm can be of type N_MAXIMUM_NORM or N_EUKLID_NORM.
567 * All arrays must have equal sizes and offsets.
568 * The complete data array inclusively offsets is used for norm calucaltion.
569 * Only non-null values are used to calculate the norm.
570 *
571 * \param a N_array_3d *
572 * \param b N_array_3d *
573 * \param type the type of the norm -> N_MAXIMUM_NORM, N_EUKLID_NORM
574 * \return double the calculated norm
575 * */
576double N_norm_array_3d(N_array_3d * a, N_array_3d * b, int type)
577{
578 int i = 0;
579 double norm = 0.0, tmp = 0.0;
580 double v1 = 0.0, v2 = 0.0;
581
582 if (a->cols_intern != b->cols_intern)
583 G_fatal_error("N_norm_array_3d: the arrays are not of equal size");
584
585 if (a->rows_intern != b->rows_intern)
586 G_fatal_error("N_norm_array_3d: the arrays are not of equal size");
587
588 if (a->depths_intern != b->depths_intern)
589 G_fatal_error("N_norm_array_3d: the arrays are not of equal size");
590
591 G_debug(3, "N_norm_array_3d: norm of a and b size %i",
593
594 for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern; i++) {
595 v1 = 0.0;
596 v2 = 0.0;
597
598 if (a->type == FCELL_TYPE) {
599 if (!Rast3d_is_null_value_num((void *)&(a->fcell_array[i]), FCELL_TYPE))
600 v1 = (double)a->fcell_array[i];
601 }
602 if (a->type == DCELL_TYPE) {
603 if (!Rast3d_is_null_value_num((void *)&(a->dcell_array[i]), DCELL_TYPE))
604 v1 = (double)a->dcell_array[i];
605 }
606 if (b->type == FCELL_TYPE) {
607 if (!Rast3d_is_null_value_num((void *)&(b->fcell_array[i]), FCELL_TYPE))
608 v2 = (double)b->fcell_array[i];
609 }
610 if (b->type == DCELL_TYPE) {
611 if (!Rast3d_is_null_value_num((void *)&(b->dcell_array[i]), DCELL_TYPE))
612 v2 = (double)b->dcell_array[i];
613 }
614
615 if (type == N_MAXIMUM_NORM) {
616 tmp = fabs(v2 - v1);
617 if ((tmp > norm))
618 norm = tmp;
619 }
620 if (type == N_EUKLID_NORM) {
621 norm += fabs(v2 - v1);
622 }
623 }
624
625 return norm;
626}
627
628/*!
629 * \brief Calculate basic statistics of the N_array_3d struct
630 *
631 * Calculates the minimum, maximum, sum and the number of
632 * non null values. The array offset can be included in the statistical calculation.
633 *
634 * \param a N_array_3d * - input array
635 * \param min double* - variable to store the computed minimum
636 * \param max double* - variable to store the computed maximum
637 * \param sum double* - variable to store the computed sum
638 * \param nonull int* - variable to store the number of non null values
639 * \param withoffset - if 1 include offset values in statistic calculation, 0 otherwise
640 * \return void
641 * */
642void N_calc_array_3d_stats(N_array_3d * a, double *min, double *max,
643 double *sum, int *nonull, int withoffset)
644{
645 int i, j, k;
646 double val;
647
648 *sum = 0.0;
649 *nonull = 0;
650
651 if (withoffset == 1) {
652
653 *min =
654 (double)N_get_array_3d_d_value(a, 0 - a->offset, 0 - a->offset,
655 0 - a->offset);
656 *max =
657 (double)N_get_array_3d_d_value(a, 0 - a->offset, 0 - a->offset,
658 0 - a->offset);
659
660 for (k = 0 - a->offset; k < a->depths + a->offset; k++) {
661 for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
662 for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
663 if (!N_is_array_3d_value_null(a, i, j, k)) {
664 val = (double)N_get_array_3d_d_value(a, i, j, k);
665 if (*min > val)
666 *min = val;
667 if (*max < val)
668 *max = val;
669 *sum += val;
670 (*nonull)++;
671 }
672 }
673 }
674 }
675 }
676 else {
677
678 *min = (double)N_get_array_3d_d_value(a, 0, 0, 0);
679 *max = (double)N_get_array_3d_d_value(a, 0, 0, 0);
680
681 for (k = 0; k < a->depths; k++) {
682 for (j = 0; j < a->rows; j++) {
683 for (i = 0; i < a->cols; i++) {
684 if (!N_is_array_3d_value_null(a, i, j, k)) {
685 val = (double)N_get_array_3d_d_value(a, i, j, k);
686 if (*min > val)
687 *min = val;
688 if (*max < val)
689 *max = val;
690 *sum += val;
691 (*nonull)++;
692 }
693 }
694 }
695 }
696 }
697
698 G_debug(3,
699 "N_calc_array_3d_stats: compute array stats, min %g, max %g, sum %g, nonull %i",
700 *min, *max, *sum, *nonull);
701
702 return;
703}
704
705/*!
706 * \brief Perform calculations with two input arrays,
707 * the result is written to a third array.
708 *
709 * All arrays must have equal sizes and offsets.
710 * The complete data array inclusively offsets is used for calucaltions.
711 * Only non-null values are used. If one array value is null,
712 * the result array value will be null too.
713 * <br><br>
714 *
715 * If a division with zero is detected, the resulting arrays
716 * value will set to null and not to NaN.
717 * <br><br>
718 *
719 * The result array is optional, if the result arrays points to NULL,
720 * a new array will be allocated with the largest arrays data type
721 * (FCELL_TYPE or DCELL_TYPE) used by the input arrays.
722 * <br><br>
723 *
724 * the calculations are of the following form:
725 *
726 * <ul>
727 * <li>result = a + b -> N_ARRAY_SUM</li>
728 * <li>result = a - b -> N_ARRAY_DIF</li>
729 * <li>result = a * b -> N_ARRAY_MUL</li>
730 * <li>result = a / b -> N_ARRAY_DIV</li>
731 * </ul>
732 *
733 * \param a N_array_3d * - first input array
734 * \param b N_array_3d * - second input array
735 * \param result N_array_3d * - the optional result array
736 * \param type - the type of calculation
737 * \return N_array_3d * - the pointer to the result array
738 * */
740 N_array_3d * result, int type)
741{
742 N_array_3d *c;
743 int i, j, k, setnull = 0;
744 double va = 0.0, vb = 0.0, vc = 0.0; /*variables used for calculation */
745
746 /*Set the pointer */
747 c = result;
748
749 /*Check the array sizes */
750 if (a->cols_intern != b->cols_intern)
751 G_fatal_error("N_math_array_3d: the arrays are not of equal size");
752 if (a->rows_intern != b->rows_intern)
753 G_fatal_error("N_math_array_3d: the arrays are not of equal size");
754 if (a->depths_intern != b->depths_intern)
755 G_fatal_error("N_math_array_3d: the arrays are not of equal size");
756 if (a->offset != b->offset)
757 G_fatal_error("N_math_array_3d: the arrays have different offsets");
758
759 G_debug(3, "N_math_array_3d: mathematical calculations, size: %i",
761
762 /*if the result array is null, allocate a new one, use the
763 * largest data type of the input arrays*/
764 if (c == NULL) {
765 if (a->type == DCELL_TYPE || b->type == DCELL_TYPE) {
766 c = N_alloc_array_3d(a->cols, a->rows, a->depths, a->offset,
767 DCELL_TYPE);
768 G_debug(3, "N_math_array_3d: array of type DCELL_TYPE created");
769 }
770 else {
771 c = N_alloc_array_3d(a->cols, a->rows, a->depths, a->offset,
772 FCELL_TYPE);
773 G_debug(3, "N_math_array_3d: array of type FCELL_TYPE created");
774 }
775 }
776 else {
777 /*Check the array sizes */
778 if (a->cols_intern != c->cols_intern)
780 ("N_math_array_3d: the arrays are not of equal size");
781 if (a->rows_intern != c->rows_intern)
783 ("N_math_array_3d: the arrays are not of equal size");
784 if (a->depths_intern != c->depths_intern)
786 ("N_math_array_3d: the arrays are not of equal size");
787 if (a->offset != c->offset)
789 ("N_math_array_3d: the arrays have different offsets");
790 }
791
792 for (k = 0 - a->offset; k < a->depths + a->offset; k++) {
793 for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
794 for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
795 if (!N_is_array_3d_value_null(a, i, j, k) &&
796 !N_is_array_3d_value_null(a, i, j, k)) {
797 /*we always calculate internally with double values */
798 va = (double)N_get_array_3d_d_value(a, i, j, k);
799 vb = (double)N_get_array_3d_d_value(b, i, j, k);
800 vc = 0;
801 setnull = 0;
802
803 switch (type) {
804 case N_ARRAY_SUM:
805 vc = va + vb;
806 break;
807 case N_ARRAY_DIF:
808 vc = va - vb;
809 break;
810 case N_ARRAY_MUL:
811 vc = va * vb;
812 break;
813 case N_ARRAY_DIV:
814 if (vb != 0)
815 vc = va / vb;
816 else
817 setnull = 1;
818 break;
819 }
820
821 if (c->type == FCELL_TYPE) {
822 if (setnull)
823 N_put_array_3d_value_null(c, i, j, k);
824 else
825 N_put_array_3d_f_value(c, i, j, k, (float)vc);
826 }
827 if (c->type == DCELL_TYPE) {
828 if (setnull)
829 N_put_array_3d_value_null(c, i, j, k);
830 else
831 N_put_array_3d_d_value(c, i, j, k, vc);
832 }
833 }
834 else {
835 N_put_array_3d_value_null(c, i, j, k);
836 }
837 }
838 }
839 }
840
841 return c;
842}
843
844/*!
845 * \brief Convert all null values to zero values
846 *
847 * The complete data array inclusively offsets is used.
848 *
849 * \param a N_array_3d *
850 * \return int - number of replaced null values
851 * */
853{
854 int i = 0, count = 0;
855
856 G_debug(3, "N_convert_array_3d_null_to_zero: convert array of size %i",
858
859 if (a->type == FCELL_TYPE)
860 for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern;
861 i++) {
862 if (Rast3d_is_null_value_num((void *)&(a->fcell_array[i]), FCELL_TYPE)) {
863 a->fcell_array[i] = 0.0;
864 count++;
865 }
866 }
867
868 if (a->type == DCELL_TYPE)
869 for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern;
870 i++) {
871 if (Rast3d_is_null_value_num((void *)&(a->dcell_array[i]), DCELL_TYPE)) {
872 a->dcell_array[i] = 0.0;
873 count++;
874 }
875 }
876
877
878 if (a->type == FCELL_TYPE)
879 G_debug(3,
880 "N_convert_array_3d_null_to_zero: %i values of type FCELL_TYPE are converted",
881 count);
882
883 if (a->type == DCELL_TYPE)
884 G_debug(3,
885 "N_convert_array_3d_null_to_zero: %i values of type DCELL_TYPE are converted",
886 count);
887
888 return count;
889}
#define N_EUKLID_NORM
Definition: N_pde.h:46
#define N_ARRAY_MUL
Definition: N_pde.h:50
#define N_MAXIMUM_NORM
Definition: N_pde.h:45
#define N_ARRAY_DIF
Definition: N_pde.h:49
#define N_ARRAY_SUM
Definition: N_pde.h:48
#define N_ARRAY_DIV
Definition: N_pde.h:51
#define NULL
Definition: ccmath.h:32
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
double b
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
int count
const char * source
Definition: lz4.h:575
void N_put_array_3d_value_null(N_array_3d *data, int col, int row, int depth)
This function writes a null value to the N_array_3d data at position col, row, depth.
Definition: n_arrays.c:1076
N_array_3d * N_alloc_array_3d(int cols, int rows, int depths, int offset, int type)
Allocate memory for a N_array_3d data structure.
Definition: n_arrays.c:726
void N_put_array_2d_f_value(N_array_2d *data, int col, int row, FCELL value)
Writes a FCELL value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:554
void N_put_array_3d_f_value(N_array_3d *data, int col, int row, int depth, float value)
This function writes a float value to the N_array_3d data at position col, row, depth.
Definition: n_arrays.c:1148
DCELL N_get_array_2d_d_value(N_array_2d *data, int col, int row)
Returns the value of type DCELL at position col, row.
Definition: n_arrays.c:378
int N_is_array_3d_value_null(N_array_3d *data, int col, int row, int depth)
This function returns 1 if value of N_array_3d data at position col, row, depth is of type null,...
Definition: n_arrays.c:881
int N_is_array_2d_value_null(N_array_2d *data, int col, int row)
Returns 1 if the value of N_array_2d struct at position col, row is of type null, otherwise 0.
Definition: n_arrays.c:231
void N_put_array_3d_d_value(N_array_3d *data, int col, int row, int depth, double value)
Writes a double value to the N_array_3d struct at position col, row, depth.
Definition: n_arrays.c:1175
N_array_2d * N_alloc_array_2d(int cols, int rows, int offset, int type)
Allocate memory for a N_array_2d data structure.
Definition: n_arrays.c:72
double N_get_array_3d_d_value(N_array_3d *data, int col, int row, int depth)
This function returns the value of type float at position col, row, depth.
Definition: n_arrays.c:990
void N_put_array_2d_value_null(N_array_2d *data, int col, int row)
Writes the null value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:458
void N_put_array_2d_c_value(N_array_2d *data, int col, int row, CELL value)
Writes a CELL value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:524
void N_put_array_2d_d_value(N_array_2d *data, int col, int row, DCELL value)
Writes a DCELL value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:584
int N_convert_array_3d_null_to_zero(N_array_3d *a)
Convert all null values to zero values.
int N_convert_array_2d_null_to_zero(N_array_2d *a)
Convert all null values to zero values.
double N_norm_array_3d(N_array_3d *a, N_array_3d *b, int type)
Calculate the norm of the two input arrays.
void N_calc_array_3d_stats(N_array_3d *a, double *min, double *max, double *sum, int *nonull, int withoffset)
Calculate basic statistics of the N_array_3d struct.
void N_calc_array_2d_stats(N_array_2d *a, double *min, double *max, double *sum, int *nonull, int withoffset)
Calculate basic statistics of the N_array_2d struct.
void N_copy_array_3d(N_array_3d *source, N_array_3d *target)
Copy the source N_array_3d struct to the target N_array_3d struct.
void N_copy_array_2d(N_array_2d *source, N_array_2d *target)
Copy the source N_array_2d struct to the target N_array_2d struct.
Definition: n_arrays_calc.c:45
double N_norm_array_2d(N_array_2d *a, N_array_2d *b, int type)
Calculate the norm of the two input arrays.
N_array_3d * N_math_array_3d(N_array_3d *a, N_array_3d *b, N_array_3d *result, int type)
Perform calculations with two input arrays, the result is written to a third array.
N_array_2d * N_math_array_2d(N_array_2d *a, N_array_2d *b, N_array_2d *result, int type)
Perform calculations with two input arrays, the result is written to a third array.
#define min(a, b)
#define max(a, b)
int type
Definition: N_pde.h:133
DCELL * dcell_array
Definition: N_pde.h:139
FCELL * fcell_array
Definition: N_pde.h:138
CELL * cell_array
Definition: N_pde.h:137
int cols
Definition: N_pde.h:134
int rows_intern
Definition: N_pde.h:135
int rows
Definition: N_pde.h:134
int offset
Definition: N_pde.h:136
int cols_intern
Definition: N_pde.h:135
int cols_intern
Definition: N_pde.h:169
int type
Definition: N_pde.h:167
int offset
Definition: N_pde.h:170
int rows
Definition: N_pde.h:168
double * dcell_array
Definition: N_pde.h:172
float * fcell_array
Definition: N_pde.h:171
int depths_intern
Definition: N_pde.h:169
int rows_intern
Definition: N_pde.h:169
int depths
Definition: N_pde.h:168
int cols
Definition: N_pde.h:168