← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to lib/gpde/N_arrays_calc.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/gpde/N_arrays_calc.c
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 managment 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
 
 
19
 
#include "grass/N_pde.h"
20
 
#include "grass/glocale.h"
21
 
#include <math.h>
22
 
 
23
 
/* ******************** 2D ARRAY FUNCTIONS *********************** */
24
 
 
25
 
/*!
26
 
 * \brief Copy the source N_array_2d struct to the target N_array_2d struct
27
 
 *
28
 
 * The arrays must have the same size and the same offset.
29
 
 *
30
 
 * The array types can be mixed, the values are automatically casted
31
 
 * and the null values are set accordingly.
32
 
 * <br><br>
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
35
 
 * <br><br>
36
 
 * This function can be called in a parallel region defined with OpenMP.
37
 
 * The copy 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
 
 * */
43
 
void N_copy_array_2d(N_array_2d * source, N_array_2d * target)
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
52
 
                ("N_copy_array_2d: the arrays are not of equal size");
53
 
 
54
 
        if (source->rows_intern != target->rows_intern)
55
 
            G_fatal_error
56
 
                ("N_copy_array_2d: the arrays are not of equal size");
57
 
 
58
 
        G_debug(3,
59
 
                "N_copy_array_2d: copy source array to target array size %i",
60
 
                source->cols_intern * source->rows_intern);
61
 
    }
62
 
 
63
 
#pragma omp for
64
 
    for (i = 0; i < source->cols_intern * source->rows_intern; i++) {
65
 
        null = 0;
66
 
        if (source->type == CELL_TYPE) {
67
 
            if (G_is_c_null_value((void *)&source->cell_array[i]))
68
 
                null = 1;
69
 
 
70
 
            if (target->type == CELL_TYPE) {
71
 
                target->cell_array[i] = source->cell_array[i];
72
 
            }
73
 
            if (target->type == FCELL_TYPE) {
74
 
                if (null)
75
 
                    G_set_f_null_value((void *)&(target->fcell_array[i]), 1);
76
 
                else
77
 
                    target->fcell_array[i] = (FCELL) source->cell_array[i];
78
 
            }
79
 
            if (target->type == DCELL_TYPE) {
80
 
                if (null)
81
 
                    G_set_d_null_value((void *)&(target->dcell_array[i]), 1);
82
 
                else
83
 
                    target->dcell_array[i] = (DCELL) source->cell_array[i];
84
 
            }
85
 
 
86
 
        }
87
 
        if (source->type == FCELL_TYPE) {
88
 
            if (G_is_f_null_value((void *)&source->fcell_array[i]))
89
 
                null = 1;
90
 
 
91
 
            if (target->type == CELL_TYPE) {
92
 
                if (null)
93
 
                    G_set_c_null_value((void *)&(target->cell_array[i]), 1);
94
 
                else
95
 
                    target->cell_array[i] = (CELL) source->fcell_array[i];
96
 
            }
97
 
            if (target->type == FCELL_TYPE) {
98
 
                target->fcell_array[i] = source->fcell_array[i];
99
 
            }
100
 
            if (target->type == DCELL_TYPE) {
101
 
                if (null)
102
 
                    G_set_d_null_value((void *)&(target->dcell_array[i]), 1);
103
 
                else
104
 
                    target->dcell_array[i] = (DCELL) source->fcell_array[i];
105
 
            }
106
 
        }
107
 
        if (source->type == DCELL_TYPE) {
108
 
            if (G_is_d_null_value((void *)&source->dcell_array[i]))
109
 
                null = 1;
110
 
 
111
 
            if (target->type == CELL_TYPE) {
112
 
                if (null)
113
 
                    G_set_c_null_value((void *)&(target->cell_array[i]), 1);
114
 
                else
115
 
                    target->cell_array[i] = (CELL) source->dcell_array[i];
116
 
            }
117
 
            if (target->type == FCELL_TYPE) {
118
 
                if (null)
119
 
                    G_set_f_null_value((void *)&(target->fcell_array[i]), 1);
120
 
                else
121
 
                    target->fcell_array[i] = (FCELL) source->dcell_array[i];
122
 
            }
123
 
            if (target->type == DCELL_TYPE) {
124
 
                target->dcell_array[i] = source->dcell_array[i];
125
 
            }
126
 
        }
127
 
    }
128
 
 
129
 
    return;
130
 
}
131
 
 
132
 
/*!
133
 
 * \brief Calculate the norm of the two input arrays
134
 
 *
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.
139
 
 *
140
 
 
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
145
 
 * */
146
 
double N_norm_array_2d(N_array_2d * a, N_array_2d * b, int type)
147
 
{
148
 
    int i = 0;
149
 
    double norm = 0.0, tmp = 0.0;
150
 
    double v1 = 0.0, v2 = 0.0;
151
 
 
152
 
    if (a->cols_intern != b->cols_intern)
153
 
        G_fatal_error("N_norm_array_2d: the arrays are not of equal size");
154
 
 
155
 
    if (a->rows_intern != b->rows_intern)
156
 
        G_fatal_error("N_norm_array_2d: the arrays are not of equal size");
157
 
 
158
 
    G_debug(3, "N_norm_array_2d: norm of a and b size %i",
159
 
            a->cols_intern * a->rows_intern);
160
 
 
161
 
    for (i = 0; i < a->cols_intern * a->rows_intern; i++) {
162
 
        v1 = 0.0;
163
 
        v2 = 0.0;
164
 
 
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];
168
 
        }
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];
172
 
        }
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];
176
 
        }
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];
180
 
        }
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];
184
 
        }
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];
188
 
        }
189
 
 
190
 
        if (type == N_MAXIMUM_NORM) {
191
 
            tmp = fabs(v2 - v1);
192
 
            if ((tmp > norm))
193
 
                norm = tmp;
194
 
        }
195
 
        if (type == N_EUKLID_NORM) {
196
 
            norm += fabs(v2 - v1);
197
 
        }
198
 
    }
199
 
 
200
 
    return norm;
201
 
}
202
 
 
203
 
/*!
204
 
 * \brief Calculate basic statistics of the N_array_2d struct 
205
 
 *
206
 
 * Calculates the minimum, maximum, sum and the number of 
207
 
 * non null values. The array offset can be included in the calculation.
208
 
 *
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 
215
 
 * \return void
216
 
 * */
217
 
void N_calc_array_2d_stats(N_array_2d * a, double *min, double *max,
218
 
                           double *sum, int *nonull, int withoffset)
219
 
{
220
 
    int i, j;
221
 
    double val;
222
 
 
223
 
    *sum = 0.0;
224
 
    *nonull = 0;
225
 
 
226
 
    if (withoffset == 1) {
227
 
 
228
 
        *min =
229
 
            (double)N_get_array_2d_d_value(a, 0 - a->offset, 0 - a->offset);
230
 
        *max =
231
 
            (double)N_get_array_2d_d_value(a, 0 - a->offset, 0 - a->offset);
232
 
 
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);
237
 
                    if (*min > val)
238
 
                        *min = val;
239
 
                    if (*max < val)
240
 
                        *max = val;
241
 
                    *sum += val;
242
 
                    (*nonull)++;
243
 
                }
244
 
            }
245
 
        }
246
 
    }
247
 
    else {
248
 
 
249
 
        *min = (double)N_get_array_2d_d_value(a, 0, 0);
250
 
        *max = (double)N_get_array_2d_d_value(a, 0, 0);
251
 
 
252
 
 
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);
257
 
                    if (*min > val)
258
 
                        *min = val;
259
 
                    if (*max < val)
260
 
                        *max = val;
261
 
                    *sum += val;
262
 
                    (*nonull)++;
263
 
                }
264
 
            }
265
 
        }
266
 
    }
267
 
 
268
 
    G_debug(3,
269
 
            "N_calc_array_2d_stats: compute array stats, min %g, max %g, sum %g, nonull %i",
270
 
            *min, *max, *sum, *nonull);
271
 
    return;
272
 
}
273
 
 
274
 
 
275
 
/*!
276
 
 * \brief Performe calculations with two input arrays, 
277
 
 * the result is written to a third array.
278
 
 *
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.
283
 
 * <br><br>
284
 
 * If a division with zero is detected, the resulting arrays 
285
 
 * value will set to null and not to NaN.
286
 
 * <br><br>
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.
290
 
 * <br><br>
291
 
 * the array computations can be of the following forms:
292
 
 *
293
 
 * <ul>
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>
298
 
 * </ul>
299
 
 *
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
305
 
 * */
306
 
N_array_2d *N_math_array_2d(N_array_2d * a, N_array_2d * b,
307
 
                            N_array_2d * result, int type)
308
 
{
309
 
    N_array_2d *c;
310
 
    int i, j, setnull = 0;
311
 
    double va = 0.0, vb = 0.0, vc = 0.0;        /*variables used for calculation */
312
 
 
313
 
    /*Set the pointer */
314
 
    c = result;
315
 
 
316
 
#pragma omp single
317
 
    {
318
 
        /*Check the array sizes */
319
 
        if (a->cols_intern != b->cols_intern)
320
 
            G_fatal_error
321
 
                ("N_math_array_2d: the arrays are not of equal size");
322
 
        if (a->rows_intern != b->rows_intern)
323
 
            G_fatal_error
324
 
                ("N_math_array_2d: the arrays are not of equal size");
325
 
        if (a->offset != b->offset)
326
 
            G_fatal_error
327
 
                ("N_math_array_2d: the arrays have different offsets");
328
 
 
329
 
        G_debug(3, "N_math_array_2d: mathematical calculations, size: %i",
330
 
                a->cols_intern * a->rows_intern);
331
 
 
332
 
        /*if the result array is null, allocate a new one, use the 
333
 
         * largest data type of the input arrays*/
334
 
        if (c == NULL) {
335
 
            if (a->type == DCELL_TYPE || b->type == DCELL_TYPE) {
336
 
                c = N_alloc_array_2d(a->cols, a->rows, a->offset, DCELL_TYPE);
337
 
                G_debug(3,
338
 
                        "N_math_array_2d: array of type DCELL_TYPE created");
339
 
            }
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);
342
 
                G_debug(3,
343
 
                        "N_math_array_2d: array of type FCELL_TYPE created");
344
 
            }
345
 
            else {
346
 
                c = N_alloc_array_2d(a->cols, a->rows, a->offset, CELL_TYPE);
347
 
                G_debug(3,
348
 
                        "N_math_array_2d: array of type CELL_TYPE created");
349
 
            }
350
 
        }
351
 
        else {
352
 
            /*Check the array sizes */
353
 
            if (a->cols_intern != c->cols_intern)
354
 
                G_fatal_error
355
 
                    ("N_math_array_2d: the arrays are not of equal size");
356
 
            if (a->rows_intern != c->rows_intern)
357
 
                G_fatal_error
358
 
                    ("N_math_array_2d: the arrays are not of equal size");
359
 
            if (a->offset != c->offset)
360
 
                G_fatal_error
361
 
                    ("N_math_array_2d: the arrays have different offsets");
362
 
        }
363
 
    }
364
 
 
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);
373
 
                vc = 0;
374
 
                setnull = 0;
375
 
 
376
 
                switch (type) {
377
 
                case N_ARRAY_SUM:
378
 
                    vc = va + vb;
379
 
                    break;
380
 
                case N_ARRAY_DIF:
381
 
                    vc = va - vb;
382
 
                    break;
383
 
                case N_ARRAY_MUL:
384
 
                    vc = va * vb;
385
 
                    break;
386
 
                case N_ARRAY_DIV:
387
 
                    if (vb != 0)
388
 
                        vc = va / vb;
389
 
                    else
390
 
                        setnull = 1;
391
 
                    break;
392
 
                }
393
 
 
394
 
                if (c->type == CELL_TYPE) {
395
 
                    if (setnull)
396
 
                        N_put_array_2d_value_null(c, i, j);
397
 
                    else
398
 
                        N_put_array_2d_c_value(c, i, j, (CELL) vc);
399
 
                }
400
 
                if (c->type == FCELL_TYPE) {
401
 
                    if (setnull)
402
 
                        N_put_array_2d_value_null(c, i, j);
403
 
                    else
404
 
                        N_put_array_2d_f_value(c, i, j, (FCELL) vc);
405
 
                }
406
 
                if (c->type == DCELL_TYPE) {
407
 
                    if (setnull)
408
 
                        N_put_array_2d_value_null(c, i, j);
409
 
                    else
410
 
                        N_put_array_2d_d_value(c, i, j, (DCELL) vc);
411
 
                }
412
 
 
413
 
            }
414
 
            else {
415
 
                N_put_array_2d_value_null(c, i, j);
416
 
            }
417
 
        }
418
 
    }
419
 
 
420
 
    return c;
421
 
}
422
 
 
423
 
/*!
424
 
 * \brief Convert all null values to zero values
425
 
 *
426
 
 * The complete data array inclusively offsets is used.
427
 
 * The array data types are automatically recognized.
428
 
 *
429
 
 * \param a N_array_2d *
430
 
 * \return int - number of replaced values
431
 
 * */
432
 
int N_convert_array_2d_null_to_zero(N_array_2d * a)
433
 
{
434
 
    int i = 0, count = 0;
435
 
 
436
 
    G_debug(3, "N_convert_array_2d_null_to_zero: convert array of size %i",
437
 
            a->cols_intern * a->rows_intern);
438
 
 
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;
443
 
                count++;
444
 
            }
445
 
        }
446
 
 
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;
451
 
                count++;
452
 
            }
453
 
        }
454
 
 
455
 
 
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;
460
 
                count++;
461
 
            }
462
 
        }
463
 
 
464
 
 
465
 
    if (a->type == CELL_TYPE)
466
 
        G_debug(2,
467
 
                "N_convert_array_2d_null_to_zero: %i values of type CELL_TYPE are converted",
468
 
                count);
469
 
    if (a->type == FCELL_TYPE)
470
 
        G_debug(2,
471
 
                "N_convert_array_2d_null_to_zero: %i valuess of type FCELL_TYPE are converted",
472
 
                count);
473
 
    if (a->type == DCELL_TYPE)
474
 
        G_debug(2,
475
 
                "N_convert_array_2d_null_to_zero: %i valuess of type DCELL_TYPE are converted",
476
 
                count);
477
 
 
478
 
    return count;
479
 
}
480
 
 
481
 
/* ******************** 3D ARRAY FUNCTIONS *********************** */
482
 
 
483
 
/*!
484
 
 * \brief Copy the source N_array_3d struct to the target N_array_3d struct
485
 
 *
486
 
 * The arrays must have the same size and the same offset.
487
 
 *
488
 
 * The array data types can be mixed, the values are automatically casted
489
 
 * and the null values are set accordingly.
490
 
 *
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
493
 
 *
494
 
 * \param source N_array_3d *
495
 
 * \param target N_array_3d *
496
 
 * \return void
497
 
 * */
498
 
void N_copy_array_3d(N_array_3d * source, N_array_3d * target)
499
 
{
500
 
    int i;
501
 
    int null;
502
 
 
503
 
    if (source->cols_intern != target->cols_intern)
504
 
        G_fatal_error("N_copy_array_3d: the arrays are not of equal size");
505
 
 
506
 
    if (source->rows_intern != target->rows_intern)
507
 
        G_fatal_error("N_copy_array_3d: the arrays are not of equal size");
508
 
 
509
 
    if (source->depths_intern != target->depths_intern)
510
 
        G_fatal_error("N_copy_array_3d: the arrays are not of equal size");
511
 
 
512
 
 
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);
516
 
 
517
 
    for (i = 0;
518
 
         i <
519
 
         source->cols_intern * source->rows_intern * source->depths_intern;
520
 
         i++) {
521
 
        null = 0;
522
 
        if (source->type == FCELL_TYPE) {
523
 
            if (G3d_isNullValueNum
524
 
                ((void *)&(source->fcell_array[i]), FCELL_TYPE))
525
 
                null = 1;
526
 
 
527
 
            if (target->type == FCELL_TYPE) {
528
 
                target->fcell_array[i] = source->fcell_array[i];
529
 
            }
530
 
            if (target->type == DCELL_TYPE) {
531
 
                if (null)
532
 
                    G3d_setNullValue((void *)&(target->dcell_array[i]), 1,
533
 
                                     DCELL_TYPE);
534
 
                else
535
 
                    target->dcell_array[i] = (double)source->fcell_array[i];
536
 
            }
537
 
 
538
 
        }
539
 
        if (source->type == DCELL_TYPE) {
540
 
            if (G3d_isNullValueNum
541
 
                ((void *)&(source->dcell_array[i]), DCELL_TYPE))
542
 
                null = 1;
543
 
 
544
 
            if (target->type == FCELL_TYPE) {
545
 
                if (null)
546
 
                    G3d_setNullValue((void *)&(target->fcell_array[i]), 1,
547
 
                                     FCELL_TYPE);
548
 
                else
549
 
                    target->fcell_array[i] = (float)source->dcell_array[i];
550
 
            }
551
 
            if (target->type == DCELL_TYPE) {
552
 
                target->dcell_array[i] = source->dcell_array[i];
553
 
            }
554
 
        }
555
 
    }
556
 
 
557
 
    return;
558
 
}
559
 
 
560
 
 
561
 
/*!
562
 
 * \brief Calculate the norm of the two input arrays
563
 
 *
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.
568
 
 *
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
573
 
 * */
574
 
double N_norm_array_3d(N_array_3d * a, N_array_3d * b, int type)
575
 
{
576
 
    int i = 0;
577
 
    double norm = 0.0, tmp = 0.0;
578
 
    double v1 = 0.0, v2 = 0.0;
579
 
 
580
 
    if (a->cols_intern != b->cols_intern)
581
 
        G_fatal_error("N_norm_array_3d: the arrays are not of equal size");
582
 
 
583
 
    if (a->rows_intern != b->rows_intern)
584
 
        G_fatal_error("N_norm_array_3d: the arrays are not of equal size");
585
 
 
586
 
    if (a->depths_intern != b->depths_intern)
587
 
        G_fatal_error("N_norm_array_3d: the arrays are not of equal size");
588
 
 
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);
591
 
 
592
 
    for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern; i++) {
593
 
        v1 = 0.0;
594
 
        v2 = 0.0;
595
 
 
596
 
        if (a->type == FCELL_TYPE) {
597
 
            if (!G3d_isNullValueNum((void *)&(a->fcell_array[i]), FCELL_TYPE))
598
 
                v1 = (double)a->fcell_array[i];
599
 
        }
600
 
        if (a->type == DCELL_TYPE) {
601
 
            if (!G3d_isNullValueNum((void *)&(a->dcell_array[i]), DCELL_TYPE))
602
 
                v1 = (double)a->dcell_array[i];
603
 
        }
604
 
        if (b->type == FCELL_TYPE) {
605
 
            if (!G3d_isNullValueNum((void *)&(b->fcell_array[i]), FCELL_TYPE))
606
 
                v2 = (double)b->fcell_array[i];
607
 
        }
608
 
        if (b->type == DCELL_TYPE) {
609
 
            if (!G3d_isNullValueNum((void *)&(b->dcell_array[i]), DCELL_TYPE))
610
 
                v2 = (double)b->dcell_array[i];
611
 
        }
612
 
 
613
 
        if (type == N_MAXIMUM_NORM) {
614
 
            tmp = fabs(v2 - v1);
615
 
            if ((tmp > norm))
616
 
                norm = tmp;
617
 
        }
618
 
        if (type == N_EUKLID_NORM) {
619
 
            norm += fabs(v2 - v1);
620
 
        }
621
 
    }
622
 
 
623
 
    return norm;
624
 
}
625
 
 
626
 
/*!
627
 
 * \brief Calculate basic statistics of the N_array_3d struct
628
 
 *
629
 
 * Calculates the minimum, maximum, sum and the number of 
630
 
 * non null values. The array offset can be included in the statistical calculation.
631
 
 *
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 
638
 
 * \return void
639
 
 * */
640
 
void N_calc_array_3d_stats(N_array_3d * a, double *min, double *max,
641
 
                           double *sum, int *nonull, int withoffset)
642
 
{
643
 
    int i, j, k;
644
 
    double val;
645
 
 
646
 
    *sum = 0.0;
647
 
    *nonull = 0;
648
 
 
649
 
    if (withoffset == 1) {
650
 
 
651
 
        *min =
652
 
            (double)N_get_array_3d_d_value(a, 0 - a->offset, 0 - a->offset,
653
 
                                           0 - a->offset);
654
 
        *max =
655
 
            (double)N_get_array_3d_d_value(a, 0 - a->offset, 0 - a->offset,
656
 
                                           0 - a->offset);
657
 
 
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);
663
 
                        if (*min > val)
664
 
                            *min = val;
665
 
                        if (*max < val)
666
 
                            *max = val;
667
 
                        *sum += val;
668
 
                        (*nonull)++;
669
 
                    }
670
 
                }
671
 
            }
672
 
        }
673
 
    }
674
 
    else {
675
 
 
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);
678
 
 
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);
684
 
                        if (*min > val)
685
 
                            *min = val;
686
 
                        if (*max < val)
687
 
                            *max = val;
688
 
                        *sum += val;
689
 
                        (*nonull)++;
690
 
                    }
691
 
                }
692
 
            }
693
 
        }
694
 
    }
695
 
 
696
 
    G_debug(3,
697
 
            "N_calc_array_3d_stats: compute array stats, min %g, max %g, sum %g, nonull %i",
698
 
            *min, *max, *sum, *nonull);
699
 
 
700
 
    return;
701
 
}
702
 
 
703
 
/*!
704
 
 * \brief Performe calculations with two input arrays, 
705
 
 * the result is written to a third array.
706
 
 *
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.
711
 
 * <br><br>
712
 
 *
713
 
 * If a division with zero is detected, the resulting arrays 
714
 
 * value will set to null and not to NaN.
715
 
 * <br><br>
716
 
 *
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.
720
 
 * <br><br>
721
 
 *
722
 
 * the calculations are of the following form:
723
 
 *
724
 
 * <ul>
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>
729
 
 * </ul>
730
 
 *
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
736
 
 * */
737
 
N_array_3d *N_math_array_3d(N_array_3d * a, N_array_3d * b,
738
 
                            N_array_3d * result, int type)
739
 
{
740
 
    N_array_3d *c;
741
 
    int i, j, k, setnull = 0;
742
 
    double va = 0.0, vb = 0.0, vc = 0.0;        /*variables used for calculation */
743
 
 
744
 
    /*Set the pointer */
745
 
    c = result;
746
 
 
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");
756
 
 
757
 
    G_debug(3, "N_math_array_3d: mathematical calculations, size: %i",
758
 
            a->cols_intern * a->rows_intern * a->depths_intern);
759
 
 
760
 
    /*if the result array is null, allocate a new one, use the 
761
 
     * largest data type of the input arrays*/
762
 
    if (c == NULL) {
763
 
        if (a->type == DCELL_TYPE || b->type == DCELL_TYPE) {
764
 
            c = N_alloc_array_3d(a->cols, a->rows, a->depths, a->offset,
765
 
                                 DCELL_TYPE);
766
 
            G_debug(3, "N_math_array_3d: array of type DCELL_TYPE created");
767
 
        }
768
 
        else {
769
 
            c = N_alloc_array_3d(a->cols, a->rows, a->depths, a->offset,
770
 
                                 FCELL_TYPE);
771
 
            G_debug(3, "N_math_array_3d: array of type FCELL_TYPE created");
772
 
        }
773
 
    }
774
 
    else {
775
 
        /*Check the array sizes */
776
 
        if (a->cols_intern != c->cols_intern)
777
 
            G_fatal_error
778
 
                ("N_math_array_3d: the arrays are not of equal size");
779
 
        if (a->rows_intern != c->rows_intern)
780
 
            G_fatal_error
781
 
                ("N_math_array_3d: the arrays are not of equal size");
782
 
        if (a->depths_intern != c->depths_intern)
783
 
            G_fatal_error
784
 
                ("N_math_array_3d: the arrays are not of equal size");
785
 
        if (a->offset != c->offset)
786
 
            G_fatal_error
787
 
                ("N_math_array_3d: the arrays have different offsets");
788
 
    }
789
 
 
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);
798
 
                    vc = 0;
799
 
                    setnull = 0;
800
 
 
801
 
                    switch (type) {
802
 
                    case N_ARRAY_SUM:
803
 
                        vc = va + vb;
804
 
                        break;
805
 
                    case N_ARRAY_DIF:
806
 
                        vc = va - vb;
807
 
                        break;
808
 
                    case N_ARRAY_MUL:
809
 
                        vc = va * vb;
810
 
                        break;
811
 
                    case N_ARRAY_DIV:
812
 
                        if (vb != 0)
813
 
                            vc = va / vb;
814
 
                        else
815
 
                            setnull = 1;
816
 
                        break;
817
 
                    }
818
 
 
819
 
                    if (c->type == FCELL_TYPE) {
820
 
                        if (setnull)
821
 
                            N_put_array_3d_value_null(c, i, j, k);
822
 
                        else
823
 
                            N_put_array_3d_f_value(c, i, j, k, (float)vc);
824
 
                    }
825
 
                    if (c->type == DCELL_TYPE) {
826
 
                        if (setnull)
827
 
                            N_put_array_3d_value_null(c, i, j, k);
828
 
                        else
829
 
                            N_put_array_3d_d_value(c, i, j, k, vc);
830
 
                    }
831
 
                }
832
 
                else {
833
 
                    N_put_array_3d_value_null(c, i, j, k);
834
 
                }
835
 
            }
836
 
        }
837
 
    }
838
 
 
839
 
    return c;
840
 
}
841
 
 
842
 
/*!
843
 
 * \brief Convert all null values to zero values
844
 
 *
845
 
 * The complete data array inclusively offsets is used.
846
 
 *
847
 
 * \param a N_array_3d *
848
 
 * \return int - number of replaced null values
849
 
 * */
850
 
int N_convert_array_3d_null_to_zero(N_array_3d * a)
851
 
{
852
 
    int i = 0, count = 0;
853
 
 
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);
856
 
 
857
 
    if (a->type == FCELL_TYPE)
858
 
        for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern;
859
 
             i++) {
860
 
            if (G3d_isNullValueNum((void *)&(a->fcell_array[i]), FCELL_TYPE)) {
861
 
                a->fcell_array[i] = 0.0;
862
 
                count++;
863
 
            }
864
 
        }
865
 
 
866
 
    if (a->type == DCELL_TYPE)
867
 
        for (i = 0; i < a->cols_intern * a->rows_intern * a->depths_intern;
868
 
             i++) {
869
 
            if (G3d_isNullValueNum((void *)&(a->dcell_array[i]), DCELL_TYPE)) {
870
 
                a->dcell_array[i] = 0.0;
871
 
                count++;
872
 
            }
873
 
        }
874
 
 
875
 
 
876
 
    if (a->type == FCELL_TYPE)
877
 
        G_debug(3,
878
 
                "N_convert_array_3d_null_to_zero: %i values of type FCELL_TYPE are converted",
879
 
                count);
880
 
 
881
 
    if (a->type == DCELL_TYPE)
882
 
        G_debug(3,
883
 
                "N_convert_array_3d_null_to_zero: %i values of type DCELL_TYPE are converted",
884
 
                count);
885
 
 
886
 
    return count;
887
 
}