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