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: 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
19 #include <math.h>
20
21 #include <grass/N_pde.h>
22 #include <grass/raster.h>
23 #include <grass/glocale.h>
24
25
26 /* ******************** 2D ARRAY FUNCTIONS *********************** */
27
28 /*!
29 * \brief Allocate memory for a N_array_2d data structure.
30 *
31 * This function allocates memory for an array of type N_array_2d
32 * and returns the pointer to the new allocated memory.
33 * <br><br>
34 * The data type of this array is set by "type" and must be
35 * CELL_TYPE, FCELL_TYPE or DCELL_TYPE accordingly to the raster map data types.
36 * The offset sets the number of boundary cols and rows.
37 * This option is useful to generate homogeneous Neumann boundary conditions around
38 * an array or to establish overlapping boundaries. The array is initialized with 0 by default.
39 * <br><br>
40 * If the offset is greater then 0, negative indices are possible.
41 * <br><br>
42 *
43 * The data structure of a array with 3 rows and cols and an offset of 1
44 * will looks like this:
45 * <br><br>
46 *
47 \verbatim
48 0 0 0 0 0
49 0 0 1 2 0
50 0 3 4 5 0
51 0 6 7 8 0
52 0 0 0 0 0
53 \endverbatim
54 *
55 * 0 is the boundary.
56 * <br><br>
57 * Internal a one dimensional array is allocated to save memory and to speed up the memory access.
58 * To access the one dimensional array with a two dimensional index use the provided
59 * get and put functions. The internal representation of the above data will look like this:
60 *
61 \verbatim
62 0 0 0 0 0 0 0 1 2 0 0 3 4 5 0 0 6 7 8 0 0 0 0 0 0
63 \endverbatim
64 *
65 * \param cols int
66 * \param rows int
67 * \param offset int
68 * \param type int
69 * \return N_array_2d *
70 *
71 * */
N_alloc_array_2d(int cols,int rows,int offset,int type)72 N_array_2d *N_alloc_array_2d(int cols, int rows, int offset, int type)
73 {
74 N_array_2d *data = NULL;
75
76 if (rows < 1 || cols < 1)
77 G_fatal_error("N_alloc_array_2d: cols and rows should be > 0");
78
79 if (type != CELL_TYPE && type != FCELL_TYPE && type != DCELL_TYPE)
80 G_fatal_error
81 ("N_alloc_array_2d: Wrong data type, should be CELL_TYPE, FCELL_TYPE or DCELL_TYPE");
82
83 data = (N_array_2d *) G_calloc(1, sizeof(N_array_2d));
84
85 data->cols = cols;
86 data->rows = rows;
87 data->type = type;
88 data->offset = offset;
89 data->rows_intern = rows + 2 * offset; /*offset position at booth sides */
90 data->cols_intern = cols + 2 * offset; /*offset position at booth sides */
91 data->cell_array = NULL;
92 data->fcell_array = NULL;
93 data->dcell_array = NULL;
94
95 if (data->type == CELL_TYPE) {
96 data->cell_array =
97 (CELL *) G_calloc((size_t) data->rows_intern * data->cols_intern,
98 sizeof(CELL));
99 G_debug(3,
100 "N_alloc_array_2d: CELL array allocated rows_intern %i cols_intern %i offset %i",
101 data->rows_intern, data->cols_intern, data->offset = offset);
102 }
103 else if (data->type == FCELL_TYPE) {
104 data->fcell_array =
105 (FCELL *) G_calloc((size_t) data->rows_intern * data->cols_intern,
106 sizeof(FCELL));
107 G_debug(3,
108 "N_alloc_array_2d: FCELL array allocated rows_intern %i cols_intern %i offset %i",
109 data->rows_intern, data->cols_intern, data->offset = offset);
110
111 }
112 else if (data->type == DCELL_TYPE) {
113 data->dcell_array =
114 (DCELL *) G_calloc((size_t) data->rows_intern * data->cols_intern,
115 sizeof(DCELL));
116 G_debug(3,
117 "N_alloc_array_2d: DCELL array allocated rows_intern %i cols_intern %i offset %i",
118 data->rows_intern, data->cols_intern, data->offset = offset);
119 }
120
121 return data;
122 }
123
124 /*!
125 * \brief Release the memory of a N_array_2d structure
126 *
127 * \param data N_array_2d *
128 * \return void
129 * */
N_free_array_2d(N_array_2d * data)130 void N_free_array_2d(N_array_2d * data)
131 {
132
133 if (data != NULL) {
134 G_debug(3, "N_free_array_2d: free N_array_2d");
135
136 if (data->type == CELL_TYPE && data->cell_array != NULL) {
137 G_free(data->cell_array);
138 }
139 else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
140 G_free(data->fcell_array);
141
142 }
143 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
144 G_free(data->dcell_array);
145 }
146
147 G_free(data);
148 data = NULL;
149
150 }
151
152 return;
153 }
154
155
156 /*!
157 * \brief Return the data type of the N_array_2d struct
158 *
159 * The data type can be CELL_TYPE, FCELL_TYPE or DCELL_TYPE accordingly to the raster map data types.
160 *
161 * \param array N_array_2d *
162 * \return type int
163 * */
N_get_array_2d_type(N_array_2d * array)164 int N_get_array_2d_type(N_array_2d * array)
165 {
166 return array->type;
167 }
168
169 /*!
170 * \brief Write the value of the N_array_2d struct at position col, row to value
171 *
172 * The value must be of the same type as the array. Otherwise you will risk data losses.
173 *
174 * \param data N_array_2d *
175 * \param col int
176 * \param row int
177 * \param value void * - this variable contains the array value at col, row position
178 * \return void
179 * */
180
N_get_array_2d_value(N_array_2d * data,int col,int row,void * value)181 void N_get_array_2d_value(N_array_2d * data, int col, int row, void *value)
182 {
183
184 if (data->offset == 0) {
185 if (data->type == CELL_TYPE && data->cell_array != NULL) {
186 *((CELL *) value) =
187 data->cell_array[row * data->cols_intern + col];
188 }
189 else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
190 *((FCELL *) value) =
191 data->fcell_array[row * data->cols_intern + col];
192 }
193 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
194 *((DCELL *) value) =
195 data->dcell_array[row * data->cols_intern + col];
196 }
197 }
198 else {
199 if (data->type == CELL_TYPE && data->cell_array != NULL) {
200 *((CELL *) value) =
201 data->cell_array[(row + data->offset) * data->cols_intern +
202 col + data->offset];
203 }
204 else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
205 *((FCELL *) value) =
206 data->fcell_array[(row + data->offset) * data->cols_intern +
207 col + data->offset];
208 }
209 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
210 *((DCELL *) value) =
211 data->dcell_array[(row + data->offset) * data->cols_intern +
212 col + data->offset];
213 }
214 }
215
216 return;
217 }
218
219 /*!
220 * \brief Returns 1 if the value of N_array_2d struct at position col, row
221 * is of type null, otherwise 0
222 *
223 * This function checks automatically the type of the array and checks for the
224 * data type null value.
225 *
226 * \param data N_array_2d *
227 * \param col int
228 * \param row int
229 * \return int - 1 = is null, 0 otherwise
230 * */
N_is_array_2d_value_null(N_array_2d * data,int col,int row)231 int N_is_array_2d_value_null(N_array_2d * data, int col, int row)
232 {
233
234 if (data->offset == 0) {
235 if (data->type == CELL_TYPE && data->cell_array != NULL) {
236 G_debug(6,
237 "N_is_array_2d_value_null: null value is of type CELL at pos [%i][%i]",
238 col, row);
239 return Rast_is_null_value((void *)
240 &(data->
241 cell_array[row * data->cols_intern +
242 col]), CELL_TYPE);
243 }
244 else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
245 G_debug(6,
246 "N_is_array_2d_value_null: null value is of type FCELL at pos [%i][%i]",
247 col, row);
248 return Rast_is_null_value((void *)
249 &(data->
250 fcell_array[row * data->cols_intern +
251 col]), FCELL_TYPE);
252 }
253 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
254 G_debug(6,
255 "N_is_array_2d_value_null: null value is of type DCELL at pos [%i][%i]",
256 col, row);
257 return Rast_is_null_value((void *)
258 &(data->
259 dcell_array[row * data->cols_intern +
260 col]), DCELL_TYPE);
261 }
262 }
263 else {
264 if (data->type == CELL_TYPE && data->cell_array != NULL) {
265 G_debug(6,
266 "N_is_array_2d_value_null: null value is of type CELL at pos [%i][%i]",
267 col, row);
268 return Rast_is_null_value((void *)
269 &(data->
270 cell_array[(row +
271 data->offset) *
272 data->cols_intern + col +
273 data->offset]), CELL_TYPE);
274 }
275 else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
276 G_debug(6,
277 "N_is_array_2d_value_null: null value is of type FCELL at pos [%i][%i]",
278 col, row);
279 return Rast_is_null_value((void *)
280 &(data->
281 fcell_array[(row +
282 data->offset) *
283 data->cols_intern + col +
284 data->offset]), FCELL_TYPE);
285 }
286 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
287 G_debug(6,
288 "N_is_array_2d_value_null: null value is of type DCELL at pos [%i][%i]",
289 col, row);
290 return Rast_is_null_value((void *)
291 &(data->
292 dcell_array[(row +
293 data->offset) *
294 data->cols_intern + col +
295 data->offset]), DCELL_TYPE);
296 }
297 }
298
299 return 0;
300 }
301
302
303 /*!
304 * \brief Returns the value of type CELL at position col, row
305 *
306 * The data array can be of type CELL, FCELL or DCELL, the value will be casted to the CELL type.
307 *
308 * \param data N_array_2d *
309 * \param col int
310 * \param row int
311 * \return CELL
312 *
313 * */
N_get_array_2d_c_value(N_array_2d * data,int col,int row)314 CELL N_get_array_2d_c_value(N_array_2d * data, int col, int row)
315 {
316 CELL value = 0;
317 FCELL fvalue = 0.0;
318 DCELL dvalue = 0.0;
319
320 switch (data->type) {
321 case CELL_TYPE:
322 N_get_array_2d_value(data, col, row, (void *)&value);
323 return (CELL) value;
324 case FCELL_TYPE:
325 N_get_array_2d_value(data, col, row, (void *)&fvalue);
326 return (CELL) fvalue;
327 case DCELL_TYPE:
328 N_get_array_2d_value(data, col, row, (void *)&dvalue);
329 return (CELL) dvalue;
330 }
331
332 return value;
333 }
334
335 /*!
336 * \brief Returns the value of type FCELL at position col, row
337 *
338 * The data array can be of type CELL, FCELL or DCELL, the value will be casted to the FCELL type.
339 *
340 * \param data N_array_2d *
341 * \param col int
342 * \param row int
343 * \return FCELL
344
345 * */
N_get_array_2d_f_value(N_array_2d * data,int col,int row)346 FCELL N_get_array_2d_f_value(N_array_2d * data, int col, int row)
347 {
348 CELL value = 0;
349 FCELL fvalue = 0.0;
350 DCELL dvalue = 0.0;
351
352 switch (data->type) {
353 case CELL_TYPE:
354 N_get_array_2d_value(data, col, row, (void *)&value);
355 return (FCELL) value;
356 case FCELL_TYPE:
357 N_get_array_2d_value(data, col, row, (void *)&fvalue);
358 return (FCELL) fvalue;
359 case DCELL_TYPE:
360 N_get_array_2d_value(data, col, row, (void *)&dvalue);
361 return (FCELL) dvalue;
362 }
363
364 return fvalue;
365 }
366
367 /*!
368 * \brief Returns the value of type DCELL at position col, row
369 *
370 * The data array can be of type CELL, FCELL or DCELL, the value will be casted to the DCELL type.
371 *
372 * \param data N_array_2d *
373 * \param col int
374 * \param row int
375 * \return DCELL
376 *
377 * */
N_get_array_2d_d_value(N_array_2d * data,int col,int row)378 DCELL N_get_array_2d_d_value(N_array_2d * data, int col, int row)
379 {
380 CELL value = 0;
381 FCELL fvalue = 0.0;
382 DCELL dvalue = 0.0;
383
384 switch (data->type) {
385 case CELL_TYPE:
386 N_get_array_2d_value(data, col, row, (void *)&value);
387 return (DCELL) value;
388 case FCELL_TYPE:
389 N_get_array_2d_value(data, col, row, (void *)&fvalue);
390 return (DCELL) fvalue;
391 case DCELL_TYPE:
392 N_get_array_2d_value(data, col, row, (void *)&dvalue);
393 return (DCELL) dvalue;
394 }
395
396 return dvalue;
397
398 }
399
400 /*!
401 * \brief Writes a value to the N_array_2d struct at position col, row
402 *
403 * The value will be automatically cast to the array type.
404 *
405 * \param data N_array_2d *
406 * \param col int
407 * \param row int
408 * \param value char *
409 * \return void
410 * */
N_put_array_2d_value(N_array_2d * data,int col,int row,char * value)411 void N_put_array_2d_value(N_array_2d * data, int col, int row, char *value)
412 {
413
414 G_debug(6, "N_put_array_2d_value: put value to array");
415
416 if (data->offset == 0) {
417 if (data->type == CELL_TYPE && data->cell_array != NULL) {
418 data->cell_array[row * data->cols_intern + col] =
419 *((CELL *) value);
420 }
421 else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
422 data->fcell_array[row * data->cols_intern + col] =
423 *((FCELL *) value);
424 }
425 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
426 data->dcell_array[row * data->cols_intern + col] =
427 *((DCELL *) value);
428 }
429 }
430 else {
431 if (data->type == CELL_TYPE && data->cell_array != NULL) {
432 data->cell_array[(row + data->offset) * data->cols_intern + col +
433 data->offset] = *((CELL *) value);
434 }
435 else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
436 data->fcell_array[(row + data->offset) * data->cols_intern + col +
437 data->offset] = *((FCELL *) value);
438 }
439 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
440 data->dcell_array[(row + data->offset) * data->cols_intern + col +
441 data->offset] = *((DCELL *) value);
442 }
443 }
444
445 return;
446 }
447
448 /*!
449 * \brief Writes the null value to the N_array_2d struct at position col, row
450 *
451 * The null value will be automatically set to the array data type (CELL, FCELL or DCELL).
452 *
453 * \param data N_array_2d *
454 * \param col int
455 * \param row int
456 * \return void
457 * */
N_put_array_2d_value_null(N_array_2d * data,int col,int row)458 void N_put_array_2d_value_null(N_array_2d * data, int col, int row)
459 {
460
461 G_debug(6,
462 "N_put_array_2d_value_null: put null value to array pos [%i][%i]",
463 col, row);
464
465 if (data->offset == 0) {
466 if (data->type == CELL_TYPE && data->cell_array != NULL) {
467 Rast_set_c_null_value((void *)
468 &(data->
469 cell_array[row * data->cols_intern + col]),
470 1);
471 }
472 else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
473 Rast_set_f_null_value((void *)
474 &(data->
475 fcell_array[row * data->cols_intern + col]),
476 1);
477 }
478 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
479 Rast_set_d_null_value((void *)
480 &(data->
481 dcell_array[row * data->cols_intern + col]),
482 1);
483 }
484 }
485 else {
486 if (data->type == CELL_TYPE && data->cell_array != NULL) {
487 Rast_set_c_null_value((void *)
488 &(data->
489 cell_array[(row +
490 data->offset) *
491 data->cols_intern + col +
492 data->offset]), 1);
493 }
494 else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
495 Rast_set_f_null_value((void *)
496 &(data->
497 fcell_array[(row +
498 data->offset) *
499 data->cols_intern + col +
500 data->offset]), 1);
501 }
502 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
503 Rast_set_d_null_value((void *)
504 &(data->
505 dcell_array[(row +
506 data->offset) *
507 data->cols_intern + col +
508 data->offset]), 1);
509 }
510 }
511
512 return;
513 }
514
515 /*!
516 * \brief Writes a CELL value to the N_array_2d struct at position col, row
517 *
518 * \param data N_array_2d *
519 * \param col int
520 * \param row int
521 * \param value CELL
522 * \return void
523 * */
N_put_array_2d_c_value(N_array_2d * data,int col,int row,CELL value)524 void N_put_array_2d_c_value(N_array_2d * data, int col, int row, CELL value)
525 {
526 FCELL fvalue;
527 DCELL dvalue;
528
529 switch (data->type) {
530 case FCELL_TYPE:
531 fvalue = (FCELL) value;
532 N_put_array_2d_value(data, col, row, (char *)&fvalue);
533 return;
534 case DCELL_TYPE:
535 dvalue = (DCELL) value;
536 N_put_array_2d_value(data, col, row, (char *)&dvalue);
537 return;
538 }
539
540 N_put_array_2d_value(data, col, row, (char *)&value);
541
542 return;
543 }
544
545 /*!
546 * \brief Writes a FCELL value to the N_array_2d struct at position col, row
547 *
548 * \param data N_array_2d *
549 * \param col int
550 * \param row int
551 * \param value FCELL
552 * \return void
553 * */
N_put_array_2d_f_value(N_array_2d * data,int col,int row,FCELL value)554 void N_put_array_2d_f_value(N_array_2d * data, int col, int row, FCELL value)
555 {
556 CELL cvalue;
557 DCELL dvalue;
558
559 switch (data->type) {
560 case CELL_TYPE:
561 cvalue = (CELL) value;
562 N_put_array_2d_value(data, col, row, (char *)&cvalue);
563 return;
564 case DCELL_TYPE:
565 dvalue = (DCELL) value;
566 N_put_array_2d_value(data, col, row, (char *)&dvalue);
567 return;
568 }
569
570 N_put_array_2d_value(data, col, row, (char *)&value);
571
572 return;
573 }
574
575 /*!
576 * \brief Writes a DCELL value to the N_array_2d struct at position col, row
577 *
578 * \param data N_array_2d *
579 * \param col int
580 * \param row int
581 * \param value DCELL
582 * \return void
583 * */
N_put_array_2d_d_value(N_array_2d * data,int col,int row,DCELL value)584 void N_put_array_2d_d_value(N_array_2d * data, int col, int row, DCELL value)
585 {
586 CELL cvalue;
587 FCELL fvalue;
588
589 switch (data->type) {
590 case CELL_TYPE:
591 cvalue = (CELL) value;
592 N_put_array_2d_value(data, col, row, (char *)&cvalue);
593 return;
594 case FCELL_TYPE:
595 fvalue = (FCELL) value;
596 N_put_array_2d_value(data, col, row, (char *)&fvalue);
597 return;
598 }
599
600 N_put_array_2d_value(data, col, row, (char *)&value);
601
602 return;
603 }
604
605 /*!
606 * \brief This function writes the data info of the array data to stdout
607 *
608 * \param data N_array_2d *
609 * \return void
610 * */
N_print_array_2d_info(N_array_2d * data)611 void N_print_array_2d_info(N_array_2d * data)
612 {
613
614 fprintf(stdout, "N_array_2d \n");
615 fprintf(stdout, "Cols %i\n", data->cols);
616 fprintf(stdout, "Rows: %i\n", data->rows);
617 fprintf(stdout, "Array type: %i\n", data->type);
618 fprintf(stdout, "Offset: %i\n", data->offset);
619 fprintf(stdout, "Internal cols: %i\n", data->cols_intern);
620 fprintf(stdout, "Internal rows: %i\n", data->rows_intern);
621 fprintf(stdout, "CELL array pointer: %p\n", data->cell_array);
622 fprintf(stdout, "FCELL array pointer: %p\n", data->fcell_array);
623 fprintf(stdout, "DCELL array pointer: %p\n", data->dcell_array);
624
625
626 return;
627 }
628
629 /*!
630 * \brief Write info and content of the N_array_2d struct to stdout
631 *
632 * Offsets are ignored
633 *
634 * \param data N_array_2d *
635 * \return void
636 * */
N_print_array_2d(N_array_2d * data)637 void N_print_array_2d(N_array_2d * data)
638 {
639 int i, j;
640
641 N_print_array_2d_info(data);
642
643 for (j = 0 - data->offset; j < data->rows + data->offset; j++) {
644 for (i = 0 - data->offset; i < data->cols + data->offset; i++) {
645 if (data->type == CELL_TYPE)
646 fprintf(stdout, "%6d ", N_get_array_2d_c_value(data, i, j));
647 else if (data->type == FCELL_TYPE)
648 fprintf(stdout, "%6.6f ", N_get_array_2d_f_value(data, i, j));
649 else if (data->type == DCELL_TYPE)
650 printf("%6.6f ", N_get_array_2d_d_value(data, i, j));
651 }
652 fprintf(stdout, "\n");
653 }
654 fprintf(stdout, "\n");
655
656 return;
657 }
658
659
660 /* ******************** 3D ARRAY FUNCTIONS *********************** */
661
662 /*!
663 * \brief Allocate memory for a N_array_3d data structure.
664 *
665 * This functions allocates an array of type N_array_3d and returns a pointer
666 * to the new allocated memory.
667 * <br><br>
668 * The data type of this array set by "type" must be
669 * FCELL_TYPE or DCELL_TYPE accordingly to the raster3d map data types.
670 * The offsets sets the number of boundary cols, rows and depths.
671 * This option is useful to generate homogeneous Neumann boundary conditions around
672 * an array or to establish overlapping boundaries. The arrays are initialized with 0 by default.
673 * <br><br>
674 * If the offset is greater then 0, negative indices are possible.
675 * The data structure of a array with 3 depths, rows and cols and an offset of 1
676 * will looks like this:
677 *
678 \verbatim
679 0 0 0 0 0
680 0 0 0 0 0
681 0 0 0 0 0
682 0 0 0 0 0
683 0 0 0 0 0
684
685 0 0 0 0 0
686 0 0 1 2 0
687 0 3 4 5 0
688 0 6 7 8 0
689 0 0 0 0 0
690
691 0 0 0 0 0
692 0 9 10 11 0
693 0 12 13 14 0
694 0 15 16 17 0
695 0 0 0 0 0
696
697 0 0 0 0 0
698 0 18 19 20 0
699 0 21 22 23 0
700 0 24 25 26 0
701 0 0 0 0 0
702
703 0 0 0 0 0
704 0 0 0 0 0
705 0 0 0 0 0
706 0 0 0 0 0
707 0 0 0 0 0
708
709 \endverbatim
710
711 The depth counts from the bottom to the top.
712
713 * <br><br>
714 * Internal a one dimensional array is allocated to speed up the memory access.
715 * To access the dimensional array with a three dimensional indexing use the provided
716 * get and put functions.
717 *
718 * \param cols int
719 * \param rows int
720 * \param depths int
721 * \param offset int
722 * \param type int
723 * \return N_array_3d *
724 *
725 * */
N_alloc_array_3d(int cols,int rows,int depths,int offset,int type)726 N_array_3d *N_alloc_array_3d(int cols, int rows, int depths, int offset,
727 int type)
728 {
729 N_array_3d *data = NULL;
730
731 if (rows < 1 || cols < 1 || depths < 1)
732 G_fatal_error
733 ("N_alloc_array_3d: depths, cols and rows should be > 0");
734
735 if (type != DCELL_TYPE && type != FCELL_TYPE)
736 G_fatal_error
737 ("N_alloc_array_3d: Wrong data type, should be FCELL_TYPE or DCELL_TYPE");
738
739 data = (N_array_3d *) G_calloc(1, sizeof(N_array_3d));
740
741 data->cols = cols;
742 data->rows = rows;
743 data->depths = depths;
744 data->type = type;
745 data->offset = offset;
746 data->rows_intern = rows + 2 * offset;
747 data->cols_intern = cols + 2 * offset;
748 data->depths_intern = depths + 2 * offset;
749 data->fcell_array = NULL;
750 data->dcell_array = NULL;
751
752 if (data->type == FCELL_TYPE) {
753 data->fcell_array =
754 (float *)G_calloc((size_t) data->depths_intern * data->rows_intern *
755 data->cols_intern, sizeof(float));
756 G_debug(3,
757 "N_alloc_array_3d: float array allocated rows_intern %i cols_intern %i depths_intern %i offset %i",
758 data->rows_intern, data->cols_intern, data->depths_intern,
759 data->offset = offset);
760 }
761 else if (data->type == DCELL_TYPE) {
762 data->dcell_array =
763 (double *)G_calloc((size_t) data->depths_intern * data->rows_intern *
764 data->cols_intern, sizeof(double));
765 G_debug(3,
766 "N_alloc_array_3d: double array allocated rows_intern %i cols_intern %i depths_intern %i offset %i",
767 data->rows_intern, data->cols_intern, data->depths_intern,
768 data->offset = offset);
769 }
770
771 return data;
772 }
773
774 /*!
775 * \brief Release the memory of a N_array_3d
776 *
777 * \param data N_array_3d *
778 * \return void
779 * */
N_free_array_3d(N_array_3d * data)780 void N_free_array_3d(N_array_3d * data)
781 {
782
783 if (data != NULL) {
784 G_debug(3, "N_free_array_3d: free N_array_3d");
785
786 if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
787 G_free(data->fcell_array);
788 }
789 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
790 G_free(data->dcell_array);
791 }
792
793 G_free(data);
794 data = NULL;
795
796 }
797
798 return;
799 }
800
801 /*!
802 * \brief Return the data type of the N_array_3d
803 *
804 * The data type can be FCELL_TYPE and DCELL_TYPE accordingly to the raster map data types.
805 *
806 * \param array N_array_3d *
807 * \return type int -- FCELL_TYPE or DCELL_TYPE
808 * */
N_get_array_3d_type(N_array_3d * array)809 int N_get_array_3d_type(N_array_3d * array)
810 {
811 return array->type;
812 }
813
814
815 /*!
816 * \brief This function writes the value of N_array_3d data at position col, row, depth
817 * to the variable value
818 *
819 * The value must be from the same type as the array. Otherwise you will risk data losses.
820 *
821 * \param data N_array_3d *
822 * \param col int
823 * \param row int
824 * \param depth int
825 * \param value void *
826 * \return void
827 * */
828 void
N_get_array_3d_value(N_array_3d * data,int col,int row,int depth,void * value)829 N_get_array_3d_value(N_array_3d * data, int col, int row, int depth,
830 void *value)
831 {
832
833 if (data->offset == 0) {
834 if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
835 *((float *)value) =
836 data->fcell_array[depth *
837 (data->rows_intern * data->cols_intern) +
838 row * data->cols_intern + col];
839 }
840 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
841 *((double *)value) =
842 data->dcell_array[depth *
843 (data->rows_intern * data->cols_intern) +
844 row * data->cols_intern + col];
845 }
846 }
847 else {
848 if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
849 *((float *)value) =
850 data->fcell_array[(depth + data->offset) *
851 (data->rows_intern * data->cols_intern) +
852 (row + data->offset) * data->cols_intern +
853 (col + data->offset)];
854
855 }
856 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
857 *((double *)value) =
858 data->dcell_array[(depth + data->offset) *
859 (data->rows_intern * data->cols_intern) +
860 (row + data->offset) * data->cols_intern +
861 (col + data->offset)];
862 }
863 }
864
865 return;
866 }
867
868 /*!
869 * \brief This function returns 1 if value of N_array_3d data at position col, row, depth
870 * is of type null, otherwise 0
871 *
872 * This function checks automatically the type of the array and checks for the
873 * data type null value.
874 *
875 * \param data N_array_3d *
876 * \param col int
877 * \param row int
878 * \param depth int
879 * \return void
880 * */
N_is_array_3d_value_null(N_array_3d * data,int col,int row,int depth)881 int N_is_array_3d_value_null(N_array_3d * data, int col, int row, int depth)
882 {
883
884 if (data->offset == 0) {
885 if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
886 G_debug(6,
887 "N_is_array_3d_value_null: null value is of type DCELL_TYPE at pos [%i][%i][%i]",
888 depth, row, col);
889 return Rast3d_is_null_value_num((void *)
890 &(data->
891 fcell_array[depth *
892 (data->rows_intern *
893 data->cols_intern) +
894 row * data->cols_intern +
895 col]), FCELL_TYPE);
896 }
897 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
898 G_debug(6,
899 "N_is_array_3d_value_null: null value is of type DCELL_TYPE at pos [%i][%i][%i]",
900 depth, row, col);
901 return Rast3d_is_null_value_num((void *)
902 &(data->
903 dcell_array[depth *
904 (data->rows_intern *
905 data->cols_intern) +
906 row * data->cols_intern +
907 col]), DCELL_TYPE);
908 }
909 }
910 else {
911 if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
912 G_debug(6,
913 "N_is_array_3d_value_null: null value is of type DCELL_TYPE at pos [%i][%i][%i]",
914 depth, row, col);
915 return Rast3d_is_null_value_num((void *)
916 &(data->
917 fcell_array[(depth +
918 data->offset) *
919 (data->rows_intern *
920 data->cols_intern) +
921 (row + data->offset)
922 * data->cols_intern +
923 (col + data->offset)]),
924 FCELL_TYPE);
925
926 }
927 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
928 G_debug(6,
929 "N_is_array_3d_value_null: null value is of type DCELL_TYPE at pos [%i][%i][%i]",
930 depth, row, col);
931 return Rast3d_is_null_value_num((void *)
932 &(data->
933 dcell_array[(depth +
934 data->offset) *
935 (data->rows_intern *
936 data->cols_intern) +
937 (row +
938 data->offset) *
939 data->cols_intern + (col +
940 data->
941 offset)]),
942 DCELL_TYPE);
943 }
944 }
945
946 return 0;
947 }
948
949 /*!
950 * \brief This function returns the value of type float at position col, row, depth
951 *
952 * The data type can be FCELL_TYPE or DCELL_TYPE accordingly to the raster map data types.
953 *
954 * \param data N_array_3d *
955 * \param col int
956 * \param row int
957 * \param depth int
958 * \return float
959 *
960 * */
N_get_array_3d_f_value(N_array_3d * data,int col,int row,int depth)961 float N_get_array_3d_f_value(N_array_3d * data, int col, int row, int depth)
962 {
963 float fvalue = 0.0;
964 double dvalue = 0.0;
965
966 switch (data->type) {
967 case FCELL_TYPE:
968 N_get_array_3d_value(data, col, row, depth, (void *)&fvalue);
969 return (float)fvalue;
970 case DCELL_TYPE:
971 N_get_array_3d_value(data, col, row, depth, (void *)&dvalue);
972 return (float)dvalue;
973 }
974
975 return fvalue;
976 }
977
978 /*!
979 * \brief This function returns the value of type float at position col, row, depth
980 *
981 * The data type can be FCELL_TYPE or DCELL_TYPE accordingly to the raster map data types.
982 *
983 * \param data N_array_3d *
984 * \param col int
985 * \param row int
986 * \param depth int
987 * \return double
988 *
989 * */
N_get_array_3d_d_value(N_array_3d * data,int col,int row,int depth)990 double N_get_array_3d_d_value(N_array_3d * data, int col, int row, int depth)
991 {
992 float fvalue = 0.0;
993 double dvalue = 0.0;
994
995 switch (data->type) {
996
997 case FCELL_TYPE:
998 N_get_array_3d_value(data, col, row, depth, (void *)&fvalue);
999 return (double)fvalue;
1000 case DCELL_TYPE:
1001 N_get_array_3d_value(data, col, row, depth, (void *)&dvalue);
1002 return (double)dvalue;
1003 }
1004
1005 return dvalue;
1006 }
1007
1008 /*!
1009 * \brief This function writes a value to the N_array_3d data at position col, row, depth
1010 *
1011 * The value will be automatically cast to the array type.
1012 *
1013 * \param data N_array_3d *
1014 * \param col int
1015 * \param row int
1016 * \param depth int
1017 * \param value cahr *
1018 * \return void
1019 * */
1020 void
N_put_array_3d_value(N_array_3d * data,int col,int row,int depth,char * value)1021 N_put_array_3d_value(N_array_3d * data, int col, int row, int depth,
1022 char *value)
1023 {
1024
1025 G_debug(6, "N_put_array_3d_value: put value to array at pos [%i][%i][%i]",
1026 depth, row, col);
1027
1028 if (data->offset == 0) {
1029 if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
1030 data->fcell_array[depth *
1031 (data->rows_intern * data->cols_intern) +
1032 row * data->cols_intern + col]
1033 = *((float *)value);
1034 }
1035 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
1036
1037 data->dcell_array[depth *
1038 (data->rows_intern * data->cols_intern) +
1039 row * data->cols_intern + col]
1040 = *((double *)value);
1041 }
1042 }
1043 else {
1044 if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
1045 data->fcell_array[(depth + data->offset) *
1046 (data->rows_intern * data->cols_intern) + (row +
1047 data->
1048 offset)
1049 * data->cols_intern + (col + data->offset)] =
1050 *((float *)value);
1051 }
1052 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
1053 data->dcell_array[(depth + data->offset) *
1054 (data->rows_intern * data->cols_intern) + (row +
1055 data->
1056 offset)
1057 * data->cols_intern + (col + data->offset)] =
1058 *((double *)value);
1059 }
1060 }
1061
1062 return;
1063 }
1064
1065 /*!
1066 * \brief This function writes a null value to the N_array_3d data at position col, row, depth
1067 *
1068 * The null value will be automatically set to the array type.
1069 *
1070 * \param data N_array_3d *
1071 * \param col int
1072 * \param row int
1073 * \param depth int
1074 * \return void
1075 * */
N_put_array_3d_value_null(N_array_3d * data,int col,int row,int depth)1076 void N_put_array_3d_value_null(N_array_3d * data, int col, int row, int depth)
1077 {
1078
1079 G_debug(6,
1080 "N_put_array_3d_value_null: put null value to array at pos [%i][%i][%i]",
1081 depth, row, col);
1082
1083 if (data->offset == 0) {
1084 if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
1085 Rast3d_set_null_value((void *)
1086 &(data->
1087 fcell_array[depth *
1088 (data->rows_intern *
1089 data->cols_intern) +
1090 row * data->cols_intern + col]), 1,
1091 FCELL_TYPE);
1092 }
1093 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
1094 Rast3d_set_null_value((void *)
1095 &(data->
1096 dcell_array[depth *
1097 (data->rows_intern *
1098 data->cols_intern) +
1099 row * data->cols_intern + col]), 1,
1100 DCELL_TYPE);
1101 }
1102 }
1103 else {
1104 if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
1105 Rast3d_set_null_value((void *)
1106 &(data->
1107 fcell_array[(depth +
1108 data->offset) *
1109 (data->rows_intern *
1110 data->cols_intern) + (row +
1111 data->
1112 offset) *
1113 data->cols_intern + (col +
1114 data->
1115 offset)]), 1,
1116 FCELL_TYPE);
1117 }
1118 else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
1119 Rast3d_set_null_value((void *)
1120 &(data->
1121 dcell_array[(depth +
1122 data->offset) *
1123 (data->rows_intern *
1124 data->cols_intern) + (row +
1125 data->
1126 offset) *
1127 data->cols_intern + (col +
1128 data->
1129 offset)]), 1,
1130 DCELL_TYPE);
1131 }
1132 }
1133
1134 return;
1135 }
1136
1137 /*!
1138 * \brief This function writes a float value to the N_array_3d data at position col, row, depth
1139 *
1140 * \param data N_array_3d *
1141 * \param col int
1142 * \param row int
1143 * \param depth int
1144 * \param value float
1145 * \return void
1146 * */
1147 void
N_put_array_3d_f_value(N_array_3d * data,int col,int row,int depth,float value)1148 N_put_array_3d_f_value(N_array_3d * data, int col, int row, int depth,
1149 float value)
1150 {
1151 double dval;
1152
1153 if (data->type == DCELL_TYPE) {
1154 dval = (double)value;
1155 N_put_array_3d_value(data, col, row, depth, (void *)&dval);
1156 }
1157 else {
1158 N_put_array_3d_value(data, col, row, depth, (void *)&value);
1159 }
1160
1161 return;
1162 }
1163
1164 /*!
1165 * \brief Writes a double value to the N_array_3d struct at position col, row, depth
1166 *
1167 * \param data N_array_3d *
1168 * \param col int
1169 * \param row int
1170 * \param depth int
1171 * \param value double
1172 * \return void
1173 * */
1174 void
N_put_array_3d_d_value(N_array_3d * data,int col,int row,int depth,double value)1175 N_put_array_3d_d_value(N_array_3d * data, int col, int row, int depth,
1176 double value)
1177 {
1178 float fval;
1179
1180 if (data->type == FCELL_TYPE) {
1181 fval = (double)value;
1182 N_put_array_3d_value(data, col, row, depth, (void *)&fval);
1183 }
1184 else {
1185 N_put_array_3d_value(data, col, row, depth, (void *)&value);
1186 }
1187
1188 return;
1189 }
1190
1191 /*!
1192 * \brief Write the info of the array to stdout
1193 *
1194 * \param data N_array_3d *
1195 * \return void
1196 * */
N_print_array_3d_info(N_array_3d * data)1197 void N_print_array_3d_info(N_array_3d * data)
1198 {
1199
1200 fprintf(stdout, "N_array_3d \n");
1201 fprintf(stdout, "Cols %i\n", data->cols);
1202 fprintf(stdout, "Rows: %i\n", data->rows);
1203 fprintf(stdout, "Depths: %i\n", data->depths);
1204 fprintf(stdout, "Array type: %i\n", data->type);
1205 fprintf(stdout, "Offset: %i\n", data->offset);
1206 fprintf(stdout, "Internal cols: %i\n", data->cols_intern);
1207 fprintf(stdout, "Internal rows: %i\n", data->rows_intern);
1208 fprintf(stdout, "Internal depths: %i\n", data->depths_intern);
1209 fprintf(stdout, "FCELL array pointer: %p\n", data->fcell_array);
1210 fprintf(stdout, "DCELL array pointer: %p\n", data->dcell_array);
1211
1212 return;
1213 }
1214
1215 /*!
1216 * \brief Write info and content of the array data to stdout
1217 *
1218 * Offsets are ignored
1219 *
1220 * \param data N_array_2d *
1221 * \return void
1222 * */
N_print_array_3d(N_array_3d * data)1223 void N_print_array_3d(N_array_3d * data)
1224 {
1225 int i, j, k;
1226
1227 N_print_array_3d_info(data);
1228
1229 for (k = 0; k < data->depths; k++) {
1230 for (j = 0; j < data->rows; j++) {
1231 for (i = 0; i < data->cols; i++) {
1232 if (data->type == FCELL_TYPE)
1233 printf("%6.6f ", N_get_array_3d_f_value(data, i, j, k));
1234 else if (data->type == DCELL_TYPE)
1235 printf("%6.6f ", N_get_array_3d_d_value(data, i, j, k));
1236 }
1237 printf("\n");
1238 }
1239 printf("\n");
1240 }
1241 printf("\n");
1242
1243 return;
1244 }
1245