2
/*****************************************************************************
4
* MODULE: Grass PDE Numerical Library
5
* AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6
* soerengebbert <at> gmx <dot> de
8
* PURPOSE: Higher level array managment functions
9
* part of the gpde library
11
* COPYRIGHT: (C) 2000 by the GRASS Development Team
13
* This program is free software under the GNU General Public
14
* License (>=v2). Read the file COPYING that comes with GRASS
17
*****************************************************************************/
19
#include "grass/N_pde.h"
20
#include "grass/glocale.h"
23
/* ******************** 2D ARRAY FUNCTIONS *********************** */
26
* \brief Copy the source N_array_2d struct to the target N_array_2d struct
28
* The arrays must have the same size and the same offset.
30
* The array types can be mixed, the values are automatically casted
31
* and the null values are set accordingly.
33
* If you copy a cell array into a dcell array, the values are casted to dcell and
34
* the null values are converted from cell-null to dcell-null
36
* This function can be called in a parallel region defined with OpenMP.
37
* The copy loop is parallelize with a openmp for pragma.
39
* \param source N_array_2d *
40
* \param target N_array_2d *
43
void N_copy_array_2d(N_array_2d * source, N_array_2d * target)
50
if (source->cols_intern != target->cols_intern)
52
("N_copy_array_2d: the arrays are not of equal size");
54
if (source->rows_intern != target->rows_intern)
56
("N_copy_array_2d: the arrays are not of equal size");
59
"N_copy_array_2d: copy source array to target array size %i",
60
source->cols_intern * source->rows_intern);
64
for (i = 0; i < source->cols_intern * source->rows_intern; i++) {
66
if (source->type == CELL_TYPE) {
67
if (G_is_c_null_value((void *)&source->cell_array[i]))
70
if (target->type == CELL_TYPE) {
71
target->cell_array[i] = source->cell_array[i];
73
if (target->type == FCELL_TYPE) {
75
G_set_f_null_value((void *)&(target->fcell_array[i]), 1);
77
target->fcell_array[i] = (FCELL) source->cell_array[i];
79
if (target->type == DCELL_TYPE) {
81
G_set_d_null_value((void *)&(target->dcell_array[i]), 1);
83
target->dcell_array[i] = (DCELL) source->cell_array[i];
87
if (source->type == FCELL_TYPE) {
88
if (G_is_f_null_value((void *)&source->fcell_array[i]))
91
if (target->type == CELL_TYPE) {
93
G_set_c_null_value((void *)&(target->cell_array[i]), 1);
95
target->cell_array[i] = (CELL) source->fcell_array[i];
97
if (target->type == FCELL_TYPE) {
98
target->fcell_array[i] = source->fcell_array[i];
100
if (target->type == DCELL_TYPE) {
102
G_set_d_null_value((void *)&(target->dcell_array[i]), 1);
104
target->dcell_array[i] = (DCELL) source->fcell_array[i];
107
if (source->type == DCELL_TYPE) {
108
if (G_is_d_null_value((void *)&source->dcell_array[i]))
111
if (target->type == CELL_TYPE) {
113
G_set_c_null_value((void *)&(target->cell_array[i]), 1);
115
target->cell_array[i] = (CELL) source->dcell_array[i];
117
if (target->type == FCELL_TYPE) {
119
G_set_f_null_value((void *)&(target->fcell_array[i]), 1);
121
target->fcell_array[i] = (FCELL) source->dcell_array[i];
123
if (target->type == DCELL_TYPE) {
124
target->dcell_array[i] = source->dcell_array[i];
133
* \brief Calculate the norm of the two input arrays
135
* The norm can be of type N_MAXIMUM_NORM or N_EUKLID_NORM.
136
* All arrays must have equal sizes and offsets.
137
* The complete data array inclusively offsets is used for norm calucaltion.
138
* Only non-null values are used to calcualte the norm.
141
* \param a N_array_2d *
142
* \param b N_array_2d *
143
* \param type the type of the norm -> N_MAXIMUM_NORM, N_EUKLID_NORM
144
* \return double the calculated norm
146
double N_norm_array_2d(N_array_2d * a, N_array_2d * b, int type)
149
double norm = 0.0, tmp = 0.0;
150
double v1 = 0.0, v2 = 0.0;
152
if (a->cols_intern != b->cols_intern)
153
G_fatal_error("N_norm_array_2d: the arrays are not of equal size");
155
if (a->rows_intern != b->rows_intern)
156
G_fatal_error("N_norm_array_2d: the arrays are not of equal size");
158
G_debug(3, "N_norm_array_2d: norm of a and b size %i",
159
a->cols_intern * a->rows_intern);
161
for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
165
if (a->type == CELL_TYPE) {
166
if (!G_is_f_null_value((void *)&(a->cell_array[i])))
167
v1 = (double)a->cell_array[i];
169
if (a->type == FCELL_TYPE) {
170
if (!G_is_f_null_value((void *)&(a->fcell_array[i])))
171
v1 = (double)a->fcell_array[i];
173
if (a->type == DCELL_TYPE) {
174
if (!G_is_f_null_value((void *)&(a->dcell_array[i])))
175
v1 = (double)a->dcell_array[i];
177
if (b->type == CELL_TYPE) {
178
if (!G_is_f_null_value((void *)&(b->cell_array[i])))
179
v2 = (double)b->cell_array[i];
181
if (b->type == FCELL_TYPE) {
182
if (!G_is_f_null_value((void *)&(b->fcell_array[i])))
183
v2 = (double)b->fcell_array[i];
185
if (b->type == DCELL_TYPE) {
186
if (!G_is_f_null_value((void *)&(b->dcell_array[i])))
187
v2 = (double)b->dcell_array[i];
190
if (type == N_MAXIMUM_NORM) {
195
if (type == N_EUKLID_NORM) {
196
norm += fabs(v2 - v1);
204
* \brief Calculate basic statistics of the N_array_2d struct
206
* Calculates the minimum, maximum, sum and the number of
207
* non null values. The array offset can be included in the calculation.
209
* \param a N_array_2d * - input array
210
* \param min double* - variable to store the computed minimum
211
* \param max double* - variable to store the computed maximum
212
* \param sum double* - variable to store the computed sum
213
* \param nonull int* - variable to store the number of non null values
214
* \param withoffset - if 1 include offset values in statistic calculation, 0 otherwise
217
void N_calc_array_2d_stats(N_array_2d * a, double *min, double *max,
218
double *sum, int *nonull, int withoffset)
226
if (withoffset == 1) {
229
(double)N_get_array_2d_d_value(a, 0 - a->offset, 0 - a->offset);
231
(double)N_get_array_2d_d_value(a, 0 - a->offset, 0 - a->offset);
233
for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
234
for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
235
if (!N_is_array_2d_value_null(a, i, j)) {
236
val = (double)N_get_array_2d_d_value(a, i, j);
249
*min = (double)N_get_array_2d_d_value(a, 0, 0);
250
*max = (double)N_get_array_2d_d_value(a, 0, 0);
253
for (j = 0; j < a->rows; j++) {
254
for (i = 0; i < a->cols; i++) {
255
if (!N_is_array_2d_value_null(a, i, j)) {
256
val = (double)N_get_array_2d_d_value(a, i, j);
269
"N_calc_array_2d_stats: compute array stats, min %g, max %g, sum %g, nonull %i",
270
*min, *max, *sum, *nonull);
276
* \brief Performe calculations with two input arrays,
277
* the result is written to a third array.
279
* All arrays must have equal sizes and offsets.
280
* The complete data array inclusively offsets is used for calucaltions.
281
* Only non-null values are computed. If one array value is null,
282
* the result array value will be null too.
284
* If a division with zero is detected, the resulting arrays
285
* value will set to null and not to NaN.
287
* The result array is optional, if the result arrays points to NULL,
288
* a new array will be allocated with the largest arrays data type
289
* (CELL, FCELL or DCELL) used by the input arrays.
291
* the array computations can be of the following forms:
294
* <li>result = a + b -> N_ARRAY_SUM</li>
295
* <li>result = a - b -> N_ARRAY_DIF</li>
296
* <li>result = a * b -> N_ARRAY_MUL</li>
297
* <li>result = a / b -> N_ARRAY_DIV</li>
300
* \param a N_array_2d * - first input array
301
* \param b N_array_2d * - second input array
302
* \param result N_array_2d * - the optional result array
303
* \param type - the type of calculation
304
* \return N_array_2d * - the pointer to the result array
306
N_array_2d *N_math_array_2d(N_array_2d * a, N_array_2d * b,
307
N_array_2d * result, int type)
310
int i, j, setnull = 0;
311
double va = 0.0, vb = 0.0, vc = 0.0; /*variables used for calculation */
318
/*Check the array sizes */
319
if (a->cols_intern != b->cols_intern)
321
("N_math_array_2d: the arrays are not of equal size");
322
if (a->rows_intern != b->rows_intern)
324
("N_math_array_2d: the arrays are not of equal size");
325
if (a->offset != b->offset)
327
("N_math_array_2d: the arrays have different offsets");
329
G_debug(3, "N_math_array_2d: mathematical calculations, size: %i",
330
a->cols_intern * a->rows_intern);
332
/*if the result array is null, allocate a new one, use the
333
* largest data type of the input arrays*/
335
if (a->type == DCELL_TYPE || b->type == DCELL_TYPE) {
336
c = N_alloc_array_2d(a->cols, a->rows, a->offset, DCELL_TYPE);
338
"N_math_array_2d: array of type DCELL_TYPE created");
340
else if (a->type == FCELL_TYPE || b->type == FCELL_TYPE) {
341
c = N_alloc_array_2d(a->cols, a->rows, a->offset, FCELL_TYPE);
343
"N_math_array_2d: array of type FCELL_TYPE created");
346
c = N_alloc_array_2d(a->cols, a->rows, a->offset, CELL_TYPE);
348
"N_math_array_2d: array of type CELL_TYPE created");
352
/*Check the array sizes */
353
if (a->cols_intern != c->cols_intern)
355
("N_math_array_2d: the arrays are not of equal size");
356
if (a->rows_intern != c->rows_intern)
358
("N_math_array_2d: the arrays are not of equal size");
359
if (a->offset != c->offset)
361
("N_math_array_2d: the arrays have different offsets");
365
#pragma omp for private(va, vb, vc, setnull)
366
for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
367
for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
368
if (!N_is_array_2d_value_null(a, i, j) &&
369
!N_is_array_2d_value_null(b, i, j)) {
370
/*we always calulate internally with double values */
371
va = (double)N_get_array_2d_d_value(a, i, j);
372
vb = (double)N_get_array_2d_d_value(b, i, j);
394
if (c->type == CELL_TYPE) {
396
N_put_array_2d_value_null(c, i, j);
398
N_put_array_2d_c_value(c, i, j, (CELL) vc);
400
if (c->type == FCELL_TYPE) {
402
N_put_array_2d_value_null(c, i, j);
404
N_put_array_2d_f_value(c, i, j, (FCELL) vc);
406
if (c->type == DCELL_TYPE) {
408
N_put_array_2d_value_null(c, i, j);
410
N_put_array_2d_d_value(c, i, j, (DCELL) vc);
415
N_put_array_2d_value_null(c, i, j);
424
* \brief Convert all null values to zero values
426
* The complete data array inclusively offsets is used.
427
* The array data types are automatically recognized.
429
* \param a N_array_2d *
430
* \return int - number of replaced values
432
int N_convert_array_2d_null_to_zero(N_array_2d * a)
434
int i = 0, count = 0;
436
G_debug(3, "N_convert_array_2d_null_to_zero: convert array of size %i",
437
a->cols_intern * a->rows_intern);
439
if (a->type == CELL_TYPE)
440
for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
441
if (G_is_c_null_value((void *)&(a->cell_array[i]))) {
442
a->cell_array[i] = 0;
447
if (a->type == FCELL_TYPE)
448
for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
449
if (G_is_f_null_value((void *)&(a->fcell_array[i]))) {
450
a->fcell_array[i] = 0.0;
456
if (a->type == DCELL_TYPE)
457
for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
458
if (G_is_d_null_value((void *)&(a->dcell_array[i]))) {
459
a->dcell_array[i] = 0.0;
465
if (a->type == CELL_TYPE)
467
"N_convert_array_2d_null_to_zero: %i values of type CELL_TYPE are converted",
469
if (a->type == FCELL_TYPE)
471
"N_convert_array_2d_null_to_zero: %i valuess of type FCELL_TYPE are converted",
473
if (a->type == DCELL_TYPE)
475
"N_convert_array_2d_null_to_zero: %i valuess of type DCELL_TYPE are converted",
481
/* ******************** 3D ARRAY FUNCTIONS *********************** */
484
* \brief Copy the source N_array_3d struct to the target N_array_3d struct
486
* The arrays must have the same size and the same offset.
488
* The array data types can be mixed, the values are automatically casted
489
* and the null values are set accordingly.
491
* If you copy a float array to a double array, the values are casted to DCELL and
492
* the null values are converted from FCELL-null to DCELL-null
494
* \param source N_array_3d *
495
* \param target N_array_3d *
498
void N_copy_array_3d(N_array_3d * source, N_array_3d * target)
503
if (source->cols_intern != target->cols_intern)
504
G_fatal_error("N_copy_array_3d: the arrays are not of equal size");
506
if (source->rows_intern != target->rows_intern)
507
G_fatal_error("N_copy_array_3d: the arrays are not of equal size");
509
if (source->depths_intern != target->depths_intern)
510
G_fatal_error("N_copy_array_3d: the arrays are not of equal size");
513
G_debug(3, "N_copy_array_3d: copy source array to target array size %i",
514
source->cols_intern * source->rows_intern *
515
source->depths_intern);
519
source->cols_intern * source->rows_intern * source->depths_intern;
522
if (source->type == FCELL_TYPE) {
523
if (G3d_isNullValueNum
524
((void *)&(source->fcell_array[i]), FCELL_TYPE))
527
if (target->type == FCELL_TYPE) {
528
target->fcell_array[i] = source->fcell_array[i];
530
if (target->type == DCELL_TYPE) {
532
G3d_setNullValue((void *)&(target->dcell_array[i]), 1,
535
target->dcell_array[i] = (double)source->fcell_array[i];
539
if (source->type == DCELL_TYPE) {
540
if (G3d_isNullValueNum
541
((void *)&(source->dcell_array[i]), DCELL_TYPE))
544
if (target->type == FCELL_TYPE) {
546
G3d_setNullValue((void *)&(target->fcell_array[i]), 1,
549
target->fcell_array[i] = (float)source->dcell_array[i];
551
if (target->type == DCELL_TYPE) {
552
target->dcell_array[i] = source->dcell_array[i];
562
* \brief Calculate the norm of the two input arrays
564
* The norm can be of type N_MAXIMUM_NORM or N_EUKLID_NORM.
565
* All arrays must have equal sizes and offsets.
566
* The complete data array inclusively offsets is used for norm calucaltion.
567
* Only non-null values are used to calcualte the norm.
569
* \param a N_array_3d *
570
* \param b N_array_3d *
571
* \param type the type of the norm -> N_MAXIMUM_NORM, N_EUKLID_NORM
572
* \return double the calculated norm
574
double N_norm_array_3d(N_array_3d * a, N_array_3d * b, int type)
577
double norm = 0.0, tmp = 0.0;
578
double v1 = 0.0, v2 = 0.0;
580
if (a->cols_intern != b->cols_intern)
581
G_fatal_error("N_norm_array_3d: the arrays are not of equal size");
583
if (a->rows_intern != b->rows_intern)
584
G_fatal_error("N_norm_array_3d: the arrays are not of equal size");
586
if (a->depths_intern != b->depths_intern)
587
G_fatal_error("N_norm_array_3d: the arrays are not of equal size");
589
G_debug(3, "N_norm_array_3d: norm of a and b size %i",
590
a->cols_intern * a->rows_intern * a->depths_intern);
592
for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern; i++) {
596
if (a->type == FCELL_TYPE) {
597
if (!G3d_isNullValueNum((void *)&(a->fcell_array[i]), FCELL_TYPE))
598
v1 = (double)a->fcell_array[i];
600
if (a->type == DCELL_TYPE) {
601
if (!G3d_isNullValueNum((void *)&(a->dcell_array[i]), DCELL_TYPE))
602
v1 = (double)a->dcell_array[i];
604
if (b->type == FCELL_TYPE) {
605
if (!G3d_isNullValueNum((void *)&(b->fcell_array[i]), FCELL_TYPE))
606
v2 = (double)b->fcell_array[i];
608
if (b->type == DCELL_TYPE) {
609
if (!G3d_isNullValueNum((void *)&(b->dcell_array[i]), DCELL_TYPE))
610
v2 = (double)b->dcell_array[i];
613
if (type == N_MAXIMUM_NORM) {
618
if (type == N_EUKLID_NORM) {
619
norm += fabs(v2 - v1);
627
* \brief Calculate basic statistics of the N_array_3d struct
629
* Calculates the minimum, maximum, sum and the number of
630
* non null values. The array offset can be included in the statistical calculation.
632
* \param a N_array_3d * - input array
633
* \param min double* - variable to store the computed minimum
634
* \param max double* - variable to store the computed maximum
635
* \param sum double* - variable to store the computed sum
636
* \param nonull int* - variable to store the number of non null values
637
* \param withoffset - if 1 include offset values in statistic calculation, 0 otherwise
640
void N_calc_array_3d_stats(N_array_3d * a, double *min, double *max,
641
double *sum, int *nonull, int withoffset)
649
if (withoffset == 1) {
652
(double)N_get_array_3d_d_value(a, 0 - a->offset, 0 - a->offset,
655
(double)N_get_array_3d_d_value(a, 0 - a->offset, 0 - a->offset,
658
for (k = 0 - a->offset; k < a->depths + a->offset; k++) {
659
for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
660
for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
661
if (!N_is_array_3d_value_null(a, i, j, k)) {
662
val = (double)N_get_array_3d_d_value(a, i, j, k);
676
*min = (double)N_get_array_3d_d_value(a, 0, 0, 0);
677
*max = (double)N_get_array_3d_d_value(a, 0, 0, 0);
679
for (k = 0; k < a->depths; k++) {
680
for (j = 0; j < a->rows; j++) {
681
for (i = 0; i < a->cols; i++) {
682
if (!N_is_array_3d_value_null(a, i, j, k)) {
683
val = (double)N_get_array_3d_d_value(a, i, j, k);
697
"N_calc_array_3d_stats: compute array stats, min %g, max %g, sum %g, nonull %i",
698
*min, *max, *sum, *nonull);
704
* \brief Performe calculations with two input arrays,
705
* the result is written to a third array.
707
* All arrays must have equal sizes and offsets.
708
* The complete data array inclusively offsets is used for calucaltions.
709
* Only non-null values are used. If one array value is null,
710
* the result array value will be null too.
713
* If a division with zero is detected, the resulting arrays
714
* value will set to null and not to NaN.
717
* The result array is optional, if the result arrays points to NULL,
718
* a new array will be allocated with the largest arrays data type
719
* (FCELL_TYPE or DCELL_TYPE) used by the input arrays.
722
* the calculations are of the following form:
725
* <li>result = a + b -> N_ARRAY_SUM</li>
726
* <li>result = a - b -> N_ARRAY_DIF</li>
727
* <li>result = a * b -> N_ARRAY_MUL</li>
728
* <li>result = a / b -> N_ARRAY_DIV</li>
731
* \param a N_array_3d * - first input array
732
* \param b N_array_3d * - second input array
733
* \param result N_array_3d * - the optional result array
734
* \param type - the type of calculation
735
* \return N_array_3d * - the pointer to the result array
737
N_array_3d *N_math_array_3d(N_array_3d * a, N_array_3d * b,
738
N_array_3d * result, int type)
741
int i, j, k, setnull = 0;
742
double va = 0.0, vb = 0.0, vc = 0.0; /*variables used for calculation */
747
/*Check the array sizes */
748
if (a->cols_intern != b->cols_intern)
749
G_fatal_error("N_math_array_3d: the arrays are not of equal size");
750
if (a->rows_intern != b->rows_intern)
751
G_fatal_error("N_math_array_3d: the arrays are not of equal size");
752
if (a->depths_intern != b->depths_intern)
753
G_fatal_error("N_math_array_3d: the arrays are not of equal size");
754
if (a->offset != b->offset)
755
G_fatal_error("N_math_array_3d: the arrays have different offsets");
757
G_debug(3, "N_math_array_3d: mathematical calculations, size: %i",
758
a->cols_intern * a->rows_intern * a->depths_intern);
760
/*if the result array is null, allocate a new one, use the
761
* largest data type of the input arrays*/
763
if (a->type == DCELL_TYPE || b->type == DCELL_TYPE) {
764
c = N_alloc_array_3d(a->cols, a->rows, a->depths, a->offset,
766
G_debug(3, "N_math_array_3d: array of type DCELL_TYPE created");
769
c = N_alloc_array_3d(a->cols, a->rows, a->depths, a->offset,
771
G_debug(3, "N_math_array_3d: array of type FCELL_TYPE created");
775
/*Check the array sizes */
776
if (a->cols_intern != c->cols_intern)
778
("N_math_array_3d: the arrays are not of equal size");
779
if (a->rows_intern != c->rows_intern)
781
("N_math_array_3d: the arrays are not of equal size");
782
if (a->depths_intern != c->depths_intern)
784
("N_math_array_3d: the arrays are not of equal size");
785
if (a->offset != c->offset)
787
("N_math_array_3d: the arrays have different offsets");
790
for (k = 0 - a->offset; k < a->depths + a->offset; k++) {
791
for (j = 0 - a->offset; j < a->rows + a->offset; j++) {
792
for (i = 0 - a->offset; i < a->cols + a->offset; i++) {
793
if (!N_is_array_3d_value_null(a, i, j, k) &&
794
!N_is_array_3d_value_null(a, i, j, k)) {
795
/*we always calulate internally with double values */
796
va = (double)N_get_array_3d_d_value(a, i, j, k);
797
vb = (double)N_get_array_3d_d_value(b, i, j, k);
819
if (c->type == FCELL_TYPE) {
821
N_put_array_3d_value_null(c, i, j, k);
823
N_put_array_3d_f_value(c, i, j, k, (float)vc);
825
if (c->type == DCELL_TYPE) {
827
N_put_array_3d_value_null(c, i, j, k);
829
N_put_array_3d_d_value(c, i, j, k, vc);
833
N_put_array_3d_value_null(c, i, j, k);
843
* \brief Convert all null values to zero values
845
* The complete data array inclusively offsets is used.
847
* \param a N_array_3d *
848
* \return int - number of replaced null values
850
int N_convert_array_3d_null_to_zero(N_array_3d * a)
852
int i = 0, count = 0;
854
G_debug(3, "N_convert_array_3d_null_to_zero: convert array of size %i",
855
a->cols_intern * a->rows_intern * a->depths_intern);
857
if (a->type == FCELL_TYPE)
858
for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern;
860
if (G3d_isNullValueNum((void *)&(a->fcell_array[i]), FCELL_TYPE)) {
861
a->fcell_array[i] = 0.0;
866
if (a->type == DCELL_TYPE)
867
for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern;
869
if (G3d_isNullValueNum((void *)&(a->dcell_array[i]), DCELL_TYPE)) {
870
a->dcell_array[i] = 0.0;
876
if (a->type == FCELL_TYPE)
878
"N_convert_array_3d_null_to_zero: %i values of type FCELL_TYPE are converted",
881
if (a->type == DCELL_TYPE)
883
"N_convert_array_3d_null_to_zero: %i values of type DCELL_TYPE are converted",