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: gradient 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 <grass/N_pde.h>
20
21 /*!
22 * \brief Allocate a N_gradient_2d structure
23 *
24 * \return N_gradient_2d *
25 *
26 * */
N_alloc_gradient_2d(void)27 N_gradient_2d *N_alloc_gradient_2d(void)
28 {
29 N_gradient_2d *grad;
30
31 grad = (N_gradient_2d *) G_calloc(1, sizeof(N_gradient_2d));
32
33 return grad;
34 }
35
36 /*!
37 * \brief Free's a N_gradient_2d structure
38 *
39 * \return void
40 *
41 * */
N_free_gradient_2d(N_gradient_2d * grad)42 void N_free_gradient_2d(N_gradient_2d * grad)
43 {
44 G_free(grad);
45 grad = NULL;
46
47 return;
48 }
49
50 /*!
51 * \brief allocate and initialize a N_gradient_2d structure
52 *
53 * \param NC double - the gradient between northern and center cell
54 * \param SC double - the gradient between southern and center cell
55 * \param WC double - the gradient between western and center cell
56 * \param EC double - the gradient between eastern and center cell
57 * \return N_gradient_2d *
58 *
59 * */
N_create_gradient_2d(double NC,double SC,double WC,double EC)60 N_gradient_2d *N_create_gradient_2d(double NC, double SC, double WC,
61 double EC)
62 {
63 N_gradient_2d *grad;
64
65 G_debug(5, "N_create_gradient_2d: create N_gradient_2d");
66
67 grad = N_alloc_gradient_2d();
68
69 grad->NC = NC;
70 grad->SC = SC;
71 grad->WC = WC;
72 grad->EC = EC;
73
74 return grad;
75 }
76
77 /*!
78 * \brief copy a N_gradient_2d structure
79 *
80 * \param source - the source N_gradient_2d struct
81 * \param target - the target N_gradient_2d struct
82 * \return int - 1 success, 0 failure while copying
83 *
84 * */
N_copy_gradient_2d(N_gradient_2d * source,N_gradient_2d * target)85 int N_copy_gradient_2d(N_gradient_2d * source, N_gradient_2d * target)
86 {
87 G_debug(5, "N_copy_gradient_2d: copy N_gradient_2d");
88
89 if (!source || !target)
90 return 0;
91
92 target->NC = source->NC;
93 target->SC = source->SC;
94 target->WC = source->WC;
95 target->EC = source->EC;
96
97 return 1;
98 }
99
100 /*!
101 * \brief Return a N_gradient_2d structure calculated from the input gradient field
102 * at position [row][col]
103 *
104 * This function returns the gradient of a cell at position [row][col] from the input gradient field.
105 * Returend is a new structure of type N_gradient_2d.
106 *
107 * \param field N_gradient_field_2d * - A two dimensional gradient field
108 * \param gradient N_gradient_2d * - the gradient structure which should be filled with data, if a NULL pointer is given, a new structure will be created
109 * \param col int
110 * \param row int
111 * \return N_gradient_2d * - the new or filled gradient structure
112 *
113 *
114 * */
N_get_gradient_2d(N_gradient_field_2d * field,N_gradient_2d * gradient,int col,int row)115 N_gradient_2d *N_get_gradient_2d(N_gradient_field_2d * field,
116 N_gradient_2d * gradient, int col, int row)
117 {
118 double NC = 0, SC = 0, WC = 0, EC = 0;
119 N_gradient_2d *grad = gradient;
120
121
122 NC = N_get_array_2d_d_value(field->y_array, col, row);
123 SC = N_get_array_2d_d_value(field->y_array, col, row + 1);
124 WC = N_get_array_2d_d_value(field->x_array, col, row);
125 EC = N_get_array_2d_d_value(field->x_array, col + 1, row);
126
127 G_debug(5,
128 "N_get_gradient_2d: calculate N_gradient_2d NC %g SC %g WC %g EC %g",
129 NC, SC, WC, EC);
130
131 /*if gradient is a NULL pointer, create a new one */
132 if (!grad) {
133 grad = N_create_gradient_2d(NC, SC, WC, EC);
134 }
135 else {
136 grad->NC = NC;
137 grad->SC = SC;
138 grad->WC = WC;
139 grad->EC = EC;
140 }
141
142 return grad;
143 }
144
145 /*!
146 * \brief Allocate a N_gradient_3d structure
147 *
148 * \return N_gradient_3d *
149 *
150 * */
N_alloc_gradient_3d(void)151 N_gradient_3d *N_alloc_gradient_3d(void)
152 {
153 N_gradient_3d *grad;
154
155 grad = (N_gradient_3d *) G_calloc(1, sizeof(N_gradient_3d));
156
157 return grad;
158 }
159
160 /*!
161 * \brief Free's a N_gradient_3d structure
162 *
163 * \return void
164 *
165 * */
N_free_gradient_3d(N_gradient_3d * grad)166 void N_free_gradient_3d(N_gradient_3d * grad)
167 {
168 G_free(grad);
169 grad = NULL;
170
171 return;
172 }
173
174
175 /*!
176 * \brief allocate and initialize a N_gradient_3d structure
177 *
178 * \param NC double - the gradient between northern and center cell
179 * \param SC double - the gradient between southern and center cell
180 * \param WC double - the gradient between western and center cell
181 * \param EC double - the gradient between eastern and center cell
182 * \param TC double - the gradient between top and center cell
183 * \param BC double - the gradient between bottom and center cell
184 * \return N_gradient_3d *
185 *
186 * */
N_create_gradient_3d(double NC,double SC,double WC,double EC,double TC,double BC)187 N_gradient_3d *N_create_gradient_3d(double NC, double SC, double WC,
188 double EC, double TC, double BC)
189 {
190 N_gradient_3d *grad;
191
192 G_debug(5, "N_create_gradient_3d: create N_gradient_3d");
193
194 grad = N_alloc_gradient_3d();
195
196 grad->NC = NC;
197 grad->SC = SC;
198 grad->WC = WC;
199 grad->EC = EC;
200 grad->TC = TC;
201 grad->BC = BC;
202
203 return grad;
204 }
205
206 /*!
207 * \brief copy a N_gradient_3d structure
208 *
209 * \param source - the source N_gradient_3d struct
210 * \param target - the target N_gradient_3d struct
211 * \return int - 1 success, 0 failure while copying
212 *
213 * */
N_copy_gradient_3d(N_gradient_3d * source,N_gradient_3d * target)214 int N_copy_gradient_3d(N_gradient_3d * source, N_gradient_3d * target)
215 {
216 G_debug(5, "N_copy_gradient_3d: copy N_gradient_3d");
217
218 if (!source || !target)
219 return 0;
220
221 target->NC = source->NC;
222 target->SC = source->SC;
223 target->WC = source->WC;
224 target->EC = source->EC;
225 target->TC = source->TC;
226 target->BC = source->BC;
227
228 return 1;
229 }
230
231
232 /*!
233 * \brief Return a N_gradient_3d structure calculated from the input gradient field
234 * at position [depth][row][col]
235 *
236 * This function returns the gradient of a 3d cell at position [depth][row][col] from the input gradient field.
237 * Returned is a new structure of type N_gradient_3d.
238 *
239 * \param field N_gradient_field_3d * - A three dimensional gradient field
240 * \param gradient N_gradient_3d * - an existing gradient structure or a NULL pointer, if a NULL pointer is providet a new structure will be returned
241 * \param col int
242 * \param row int
243 * \param depth int
244 * \return N_gradient_3d *
245 *
246 *
247 * */
N_get_gradient_3d(N_gradient_field_3d * field,N_gradient_3d * gradient,int col,int row,int depth)248 N_gradient_3d *N_get_gradient_3d(N_gradient_field_3d * field,
249 N_gradient_3d * gradient, int col, int row,
250 int depth)
251 {
252 double NC, SC, WC, EC, TC, BC;
253 N_gradient_3d *grad = gradient;
254
255 NC = N_get_array_3d_d_value(field->y_array, col, row, depth);
256 SC = N_get_array_3d_d_value(field->y_array, col, row + 1, depth);
257 WC = N_get_array_3d_d_value(field->x_array, col, row, depth);
258 EC = N_get_array_3d_d_value(field->x_array, col + 1, row, depth);
259 BC = N_get_array_3d_d_value(field->z_array, col, row, depth);
260 TC = N_get_array_3d_d_value(field->z_array, col, row, depth + 1);
261
262 G_debug(6,
263 "N_get_gradient_3d: calculate N_gradient_3d NC %g SC %g WC %g EC %g TC %g BC %g",
264 NC, SC, WC, EC, TC, BC);
265
266 /*if gradient is a NULL pointer, create a new one */
267 if (!grad) {
268 grad = N_create_gradient_3d(NC, SC, WC, EC, TC, BC);
269 }
270 else {
271 grad->NC = NC;
272 grad->SC = SC;
273 grad->WC = WC;
274 grad->EC = EC;
275 grad->BC = BC;
276 grad->TC = TC;
277 }
278
279 return grad;
280 }
281
282 /*!
283 * \brief Allocate a N_gradient_neighbours_x structure
284 *
285 * This structure contains all neighbour gradients in x direction of one cell
286 *
287 * \return N_gradient_neighbours_x *
288 *
289 * */
N_alloc_gradient_neighbours_x(void)290 N_gradient_neighbours_x *N_alloc_gradient_neighbours_x(void)
291 {
292 N_gradient_neighbours_x *grad;
293
294 grad =
295 (N_gradient_neighbours_x *) G_calloc(1,
296 sizeof(N_gradient_neighbours_x));
297
298 return grad;
299 }
300
301 /*!
302 * \brief Free's a N_gradient_neighbours_x structure
303 *
304 * \return void
305 *
306 * */
N_free_gradient_neighbours_x(N_gradient_neighbours_x * grad)307 void N_free_gradient_neighbours_x(N_gradient_neighbours_x * grad)
308 {
309 G_free(grad);
310 grad = NULL;
311
312 return;
313 }
314
315
316 /*!
317 * \brief Allocate and initialize a N_gradient_neighbours_x structure
318 *
319 * \param NWN double - the gradient between north-west and northern cell
320 * \param NEN double - the gradient between north-east and northern cell
321 * \param WC double - the gradient between western and center cell
322 * \param EC double - the gradient between eastern and center cell
323 * \param SWS double - the gradient between south-west and southern cell
324 * \param SES double - the gradient between south-east and southern cell
325 * \return N_gradient_neighbours_x *
326
327 *
328 * */
N_create_gradient_neighbours_x(double NWN,double NEN,double WC,double EC,double SWS,double SES)329 N_gradient_neighbours_x *N_create_gradient_neighbours_x(double NWN,
330 double NEN, double WC,
331 double EC, double SWS,
332 double SES)
333 {
334 N_gradient_neighbours_x *grad;
335
336 G_debug(6,
337 "N_create_gradient_neighbours_x: create N_gradient_neighbours_x");
338
339 grad = N_alloc_gradient_neighbours_x();
340
341 grad->NWN = NWN;
342 grad->NEN = NEN;
343 grad->WC = WC;
344 grad->EC = EC;
345 grad->SWS = SWS;
346 grad->SES = SES;
347
348 return grad;
349 }
350
351 /*!
352 * \brief copy a N_gradient_neighbours_x structure
353 *
354 * \param source - the source N_gradient_neighbours_x struct
355 * \param target - the target N_gradient_neighbours_x struct
356 * \return int - 1 success, 0 failure while copying
357 *
358 * */
359 int
N_copy_gradient_neighbours_x(N_gradient_neighbours_x * source,N_gradient_neighbours_x * target)360 N_copy_gradient_neighbours_x(N_gradient_neighbours_x * source,
361 N_gradient_neighbours_x * target)
362 {
363 G_debug(6, "N_copy_gradient_neighbours_x: copy N_gradient_neighbours_x");
364
365 if (!source || !target)
366 return 0;
367
368 target->NWN = source->NWN;
369 target->NEN = source->NEN;
370 target->WC = source->WC;
371 target->EC = source->EC;
372 target->SWS = source->SWS;
373 target->SES = source->SES;
374
375 return 1;
376 }
377
378 /*!
379 * \brief Allocate a N_gradient_neighbours_y structure
380 *
381 * This structure contains all neighbour gradients in y direction of one cell
382 *
383 * \return N_gradient_neighbours_y *
384 *
385 * */
N_alloc_gradient_neighbours_y(void)386 N_gradient_neighbours_y *N_alloc_gradient_neighbours_y(void)
387 {
388 N_gradient_neighbours_y *grad;
389
390 grad =
391 (N_gradient_neighbours_y *) G_calloc(1,
392 sizeof(N_gradient_neighbours_y));
393
394 return grad;
395 }
396
397 /*!
398 * \brief Free's a N_gradient_neighbours_y structure
399 *
400 * \return void
401 *
402 * */
N_free_gradient_neighbours_y(N_gradient_neighbours_y * grad)403 void N_free_gradient_neighbours_y(N_gradient_neighbours_y * grad)
404 {
405 G_free(grad);
406 grad = NULL;
407
408 return;
409 }
410
411 /*!
412 * \brief Allocate and initialize a N_gradient_neighbours_y structure
413 *
414 * \param NWW double - the gradient between north-west and western cell
415 * \param NEE double - the gradient between north-east and eastern cell
416 * \param NC double - the gradient between northern and center cell
417 * \param SC double - the gradient between southern and center cell
418 * \param SWW double - the gradient between south-west and western cell
419 * \param SEE double - the gradient between south-east and eastern cell
420 * \return N_gradient_neighbours_y *
421
422 *
423 * */
N_create_gradient_neighbours_y(double NWW,double NEE,double NC,double SC,double SWW,double SEE)424 N_gradient_neighbours_y *N_create_gradient_neighbours_y(double NWW,
425 double NEE, double NC,
426 double SC, double SWW,
427 double SEE)
428 {
429 N_gradient_neighbours_y *grad;
430
431 G_debug(6,
432 "N_create_gradient_neighbours_y: create N_gradient_neighbours_y");
433
434 grad = N_alloc_gradient_neighbours_y();
435
436 grad->NWW = NWW;
437 grad->NEE = NEE;
438 grad->NC = NC;
439 grad->SC = SC;
440 grad->SWW = SWW;
441 grad->SEE = SEE;
442
443 return grad;
444 }
445
446 /*!
447 * \brief copy a N_gradient_neighbours_y structure
448 *
449 * \param source - the source N_gradient_neighbours_y struct
450 * \param target - the target N_gradient_neighbours_y struct
451 * \return int - 1 success, 0 failure while copying
452 *
453 * */
454 int
N_copy_gradient_neighbours_y(N_gradient_neighbours_y * source,N_gradient_neighbours_y * target)455 N_copy_gradient_neighbours_y(N_gradient_neighbours_y * source,
456 N_gradient_neighbours_y * target)
457 {
458 G_debug(6, "N_copy_gradient_neighbours_y: copy N_gradient_neighbours_y");
459
460 if (!source || !target)
461 return 0;
462
463 target->NWW = source->NWW;
464 target->NEE = source->NEE;
465 target->NC = source->NC;
466 target->SC = source->SC;
467 target->SWW = source->SWW;
468 target->SEE = source->SEE;
469
470 return 1;
471 }
472
473 /*!
474 * \brief Allocate a N_gradient_neighbours_z structure
475 *
476 * This structure contains all neighbour gradients in z direction of one cell
477 *
478 * \return N_gradient_neighbours_z *
479 *
480 * */
N_alloc_gradient_neighbours_z(void)481 N_gradient_neighbours_z *N_alloc_gradient_neighbours_z(void)
482 {
483 N_gradient_neighbours_z *grad;
484
485 grad =
486 (N_gradient_neighbours_z *) G_calloc(1,
487 sizeof(N_gradient_neighbours_z));
488
489 return grad;
490 }
491
492 /*!
493 * \brief Free's a N_gradient_neighbours_z structure
494 *
495 * \return void
496 *
497 * */
N_free_gradient_neighbours_z(N_gradient_neighbours_z * grad)498 void N_free_gradient_neighbours_z(N_gradient_neighbours_z * grad)
499 {
500 G_free(grad);
501 grad = NULL;
502
503 return;
504 }
505
506 /*!
507 * \brief Allocate and initialize a N_gradient_neighbours_z structure
508 *
509 * \param NWZ double - the gradient between upper and lower north-western cells
510 * \param NZ double - the gradient between upper and lower northern cells
511 * \param NEZ double - the gradient between upper and lower north-eastern cells
512 * \param WZ double - the gradient between upper and lower western cells
513 * \param CZ double - the gradient between upper and lower center cells
514 * \param EZ double - the gradient between upper and lower eastern cells
515 * \param SWZ double - the gradient between upper and lower south-western cells
516 * \param SZ double - the gradient between upper and lower southern cells
517 * \param SEZ double - the gradient between upper and lower south-eastern cells
518 * \return N_gradient_neighbours_z *
519
520 *
521 * */
N_create_gradient_neighbours_z(double NWZ,double NZ,double NEZ,double WZ,double CZ,double EZ,double SWZ,double SZ,double SEZ)522 N_gradient_neighbours_z *N_create_gradient_neighbours_z(double NWZ, double NZ,
523 double NEZ, double WZ,
524 double CZ, double EZ,
525 double SWZ, double SZ,
526 double SEZ)
527 {
528 N_gradient_neighbours_z *grad;
529
530 G_debug(6,
531 "N_create_gradient_neighbours_z: create N_gradient_neighbours_z");
532
533 grad = N_alloc_gradient_neighbours_z();
534
535 grad->NWZ = NWZ;
536 grad->NZ = NZ;
537 grad->NEZ = NEZ;
538 grad->WZ = WZ;
539 grad->CZ = CZ;
540 grad->EZ = EZ;
541 grad->SWZ = SWZ;
542 grad->SZ = SZ;
543 grad->SEZ = SEZ;
544
545 return grad;
546 }
547
548 /*!
549 * \brief copy a N_gradient_neighbours_z structure
550 *
551 * \param source - the source N_gradient_neighbours_z struct
552 * \param target - the target N_gradient_neighbours_z struct
553 * \return int - 1 success, 0 failure while copying
554 *
555 * */
556 int
N_copy_gradient_neighbours_z(N_gradient_neighbours_z * source,N_gradient_neighbours_z * target)557 N_copy_gradient_neighbours_z(N_gradient_neighbours_z * source,
558 N_gradient_neighbours_z * target)
559 {
560 G_debug(6, "N_copy_gradient_neighbours_z: copy N_gradient_neighbours_z");
561
562 if (!source || !target)
563 return 0;
564
565 target->NWZ = source->NWZ;
566 target->NZ = source->NZ;
567 target->NEZ = source->NEZ;
568 target->WZ = source->WZ;
569 target->CZ = source->CZ;
570 target->EZ = source->EZ;
571 target->SWZ = source->SWZ;
572 target->SZ = source->SZ;
573 target->SEZ = source->SEZ;
574
575 return 1;
576 }
577
578 /*!
579 * \brief Allocate a N_gradient_neighbours_2d structure
580 *
581 * This structure contains all neighbour gradients in all directions of one cell
582 * in a 2d raster layer
583 *
584 * \return N_gradient_neighbours_2d *
585 *
586 * */
N_alloc_gradient_neighbours_2d(void)587 N_gradient_neighbours_2d *N_alloc_gradient_neighbours_2d(void)
588 {
589 N_gradient_neighbours_2d *grad;
590
591 grad =
592 (N_gradient_neighbours_2d *) G_calloc(1,
593 sizeof
594 (N_gradient_neighbours_2d));
595
596 grad->x = N_alloc_gradient_neighbours_x();
597 grad->y = N_alloc_gradient_neighbours_y();
598
599 return grad;
600 }
601
602 /*!
603 * \brief Free's a N_gradient_neighbours_2d structure
604 *
605 * \return void
606 *
607 * */
N_free_gradient_neighbours_2d(N_gradient_neighbours_2d * grad)608 void N_free_gradient_neighbours_2d(N_gradient_neighbours_2d * grad)
609 {
610
611 N_free_gradient_neighbours_x(grad->x);
612 N_free_gradient_neighbours_y(grad->y);
613
614 G_free(grad);
615 grad = NULL;
616
617 return;
618 }
619
620 /*!
621 * \brief Allocate and initialize a N_gradient_neighbours_2d structure
622 *
623 * The parameter N_gradient_neighbours x and y are copied into the new allocated structure
624 * and can be deleted after the initializing
625 *
626 * \return N_gradient_neighbours_2d * -- if failure NULL is returned
627 *
628 * */
629 N_gradient_neighbours_2d
N_create_gradient_neighbours_2d(N_gradient_neighbours_x * x,N_gradient_neighbours_y * y)630 * N_create_gradient_neighbours_2d(N_gradient_neighbours_x * x,
631 N_gradient_neighbours_y * y)
632 {
633 N_gradient_neighbours_2d *grad;
634 int fail = 0;
635
636 G_debug(5,
637 "N_create_gradient_neighbours_2d: create N_gradient_neighbours_2d");
638
639 grad = N_alloc_gradient_neighbours_2d();
640
641 if (!N_copy_gradient_neighbours_x(x, grad->x))
642 fail++;
643 if (!N_copy_gradient_neighbours_y(y, grad->y))
644 fail++;
645
646 if (fail > 0) {
647 N_free_gradient_neighbours_2d(grad);
648 grad = NULL;
649 }
650
651 return grad;
652 }
653
654 /*!
655 * \brief copy a N_gradient_neighbours_2d structure
656 *
657 * \param source - the source N_gradient_neighbours_2d struct
658 * \param target - the target N_gradient_neighbours_2d struct
659 * \return int - 1 success, 0 failure while copying
660 *
661 * */
662 int
N_copy_gradient_neighbours_2d(N_gradient_neighbours_2d * source,N_gradient_neighbours_2d * target)663 N_copy_gradient_neighbours_2d(N_gradient_neighbours_2d * source,
664 N_gradient_neighbours_2d * target)
665 {
666 int fail = 0;
667
668 G_debug(5,
669 "N_copy_gradient_neighbours_2d: copy N_gradient_neighbours_2d");
670
671 if (!source || !target)
672 return 0;
673
674 if (!(N_copy_gradient_neighbours_x(source->x, target->x)))
675 fail++;
676 if (!(N_copy_gradient_neighbours_y(source->y, target->y)))
677 fail++;
678
679 if (fail > 0) {
680 return 0;
681 }
682
683 return 1;
684 }
685
686 /*!
687 * \brief Return a N_gradient_neighbours_2d structure calculated from the input gradient field
688 * at position [row][col]
689 *
690 * This function returns the gradient neighbours in x and y dierection
691 * of a cell at position [row][col] from the input gradient field.
692 * Returend is a pointer to a structure of type N_gradient_neighbours_2d.
693 *
694 * \param field N_gradient_field_2d * - A two dimensional gradient field
695 * \param gradient N_gradient_neighbours_2d * - the gradient structure which should be filled with data, if a NULL pointer is given, a new structure will be created
696 * \param col int
697 * \param row int
698 * \return N_gradient_neighbours_2d * - the new or filled gradient structure
699 *
700 *
701 * */
N_get_gradient_neighbours_2d(N_gradient_field_2d * field,N_gradient_neighbours_2d * gradient,int col,int row)702 N_gradient_neighbours_2d *N_get_gradient_neighbours_2d(N_gradient_field_2d *
703 field,
704 N_gradient_neighbours_2d
705 * gradient, int col,
706 int row)
707 {
708 double NWN, NEN, WC, EC, SWS, SES;
709 double NWW, NEE, NC, SC, SWW, SEE;
710 N_gradient_neighbours_2d *grad = NULL;
711 N_gradient_neighbours_x *grad_x = NULL;
712 N_gradient_neighbours_y *grad_y = NULL;
713
714
715 NWN = N_get_array_2d_d_value(field->x_array, col, row - 1);
716 NEN = N_get_array_2d_d_value(field->x_array, col + 1, row - 1);
717 WC = N_get_array_2d_d_value(field->x_array, col, row);
718 EC = N_get_array_2d_d_value(field->x_array, col + 1, row);
719 SWS = N_get_array_2d_d_value(field->x_array, col, row + 1);
720 SES = N_get_array_2d_d_value(field->x_array, col + 1, row + 1);
721
722 NWW = N_get_array_2d_d_value(field->y_array, col - 1, row);
723 NEE = N_get_array_2d_d_value(field->y_array, col + 1, row);
724 NC = N_get_array_2d_d_value(field->y_array, col, row);
725 SC = N_get_array_2d_d_value(field->y_array, col, row + 1);
726 SWW = N_get_array_2d_d_value(field->y_array, col - 1, row + 1);
727 SEE = N_get_array_2d_d_value(field->y_array, col + 1, row + 1);
728
729
730 grad_x = N_create_gradient_neighbours_x(NWN, NEN, WC, EC, SWS, SES);
731 grad_y = N_create_gradient_neighbours_y(NWW, NEE, NC, SC, SWW, SEE);
732
733 G_debug(5,
734 "N_get_gradient_neighbours_2d: calculate N_gradient_neighbours_x NWN %g NEN %g WC %g EC %g SWS %g SES %g",
735 NWN, NEN, WC, EC, SWS, SES);
736
737 G_debug(5,
738 "N_get_gradient_neighbours_2d: calculate N_gradient_neighbours_y NWW %g NEE %g NC %g SC %g SWW %g SEE %g",
739 NWW, NEE, NC, SC, SWW, SEE);
740
741
742 /*if gradient is a NULL pointer, create a new one */
743 if (!gradient) {
744 grad = N_create_gradient_neighbours_2d(grad_x, grad_y);
745 gradient = grad;
746 }
747 else {
748 grad = N_create_gradient_neighbours_2d(grad_x, grad_y);
749 N_copy_gradient_neighbours_2d(grad, gradient);
750 N_free_gradient_neighbours_2d(grad);
751 }
752
753 N_free_gradient_neighbours_x(grad_x);
754 N_free_gradient_neighbours_y(grad_y);
755
756 return gradient;
757 }
758
759
760 /*!
761 * \brief Allocate a N_gradient_neighbours_3d structure
762 *
763 * This structure contains all neighbour gradients in all directions of one cell
764 * in a 3d raster layer
765 *
766 * \return N_gradient_neighbours_3d *
767 *
768 * */
N_alloc_gradient_neighbours_3d(void)769 N_gradient_neighbours_3d *N_alloc_gradient_neighbours_3d(void)
770 {
771 N_gradient_neighbours_3d *grad;
772
773 grad =
774 (N_gradient_neighbours_3d *) G_calloc(1,
775 sizeof
776 (N_gradient_neighbours_3d));
777
778 grad->xt = N_alloc_gradient_neighbours_x();
779 grad->xc = N_alloc_gradient_neighbours_x();
780 grad->xb = N_alloc_gradient_neighbours_x();
781 grad->yt = N_alloc_gradient_neighbours_y();
782 grad->yc = N_alloc_gradient_neighbours_y();
783 grad->yb = N_alloc_gradient_neighbours_y();
784 grad->zt = N_alloc_gradient_neighbours_z();
785 grad->zb = N_alloc_gradient_neighbours_z();
786
787 return grad;
788 }
789
790 /*!
791 * \brief Free's a N_gradient_neighbours_3d structure
792 *
793 * \return void
794 *
795 * */
N_free_gradient_neighbours_3d(N_gradient_neighbours_3d * grad)796 void N_free_gradient_neighbours_3d(N_gradient_neighbours_3d * grad)
797 {
798
799 N_free_gradient_neighbours_x(grad->xt);
800 N_free_gradient_neighbours_x(grad->xc);
801 N_free_gradient_neighbours_x(grad->xb);
802 N_free_gradient_neighbours_y(grad->yt);
803 N_free_gradient_neighbours_y(grad->yc);
804 N_free_gradient_neighbours_y(grad->yb);
805 N_free_gradient_neighbours_z(grad->zt);
806 N_free_gradient_neighbours_z(grad->zb);
807
808 G_free(grad);
809 grad = NULL;
810
811 return;
812 }
813
814 /*!
815 * \brief Allocate and initialize a N_gradient_neighbours_3d structure
816 *
817 * The parameter N_gradient_neighbours x(tcb) and y(tcb) and z(tb) are copied into the new allocated structure
818 * and can be deleted after the initializing
819 *
820 * \return N_gradient_neighbours_3d * -- if failure NULL is returned
821
822 *
823 * */
824 N_gradient_neighbours_3d
N_create_gradient_neighbours_3d(N_gradient_neighbours_x * xt,N_gradient_neighbours_x * xc,N_gradient_neighbours_x * xb,N_gradient_neighbours_y * yt,N_gradient_neighbours_y * yc,N_gradient_neighbours_y * yb,N_gradient_neighbours_z * zt,N_gradient_neighbours_z * zb)825 * N_create_gradient_neighbours_3d(N_gradient_neighbours_x * xt,
826 N_gradient_neighbours_x * xc,
827 N_gradient_neighbours_x * xb,
828 N_gradient_neighbours_y * yt,
829 N_gradient_neighbours_y * yc,
830 N_gradient_neighbours_y * yb,
831 N_gradient_neighbours_z * zt,
832 N_gradient_neighbours_z * zb)
833 {
834 N_gradient_neighbours_3d *grad;
835 int fail = 0;
836
837 G_debug(5,
838 "N_create_gradient_neighbours_3d: create N_gradient_neighbours_3d");
839
840 grad = N_alloc_gradient_neighbours_3d();
841
842 if (!(N_copy_gradient_neighbours_x(xt, grad->xt)))
843 fail++;
844 if (!(N_copy_gradient_neighbours_x(xc, grad->xc)))
845 fail++;
846 if (!(N_copy_gradient_neighbours_x(xb, grad->xb)))
847 fail++;
848 if (!(N_copy_gradient_neighbours_y(yt, grad->yt)))
849 fail++;
850 if (!(N_copy_gradient_neighbours_y(yc, grad->yc)))
851 fail++;
852 if (!(N_copy_gradient_neighbours_y(yb, grad->yb)))
853 fail++;
854 if (!(N_copy_gradient_neighbours_z(zt, grad->zt)))
855 fail++;
856 if (!(N_copy_gradient_neighbours_z(zb, grad->zb)))
857 fail++;
858
859 if (fail > 0) {
860 return NULL;
861 }
862
863 return grad;
864 }
865
866 /*!
867 * \brief copy a N_gradient_neighbours_3d structure
868 *
869 * \param source - the source N_gradient_neighbours_3d struct
870 * \param target - the target N_gradient_neighbours_3d struct
871 * \return int - 1 success, 0 failure while copying
872 *
873 * */
874 int
N_copy_gradient_neighbours_3d(N_gradient_neighbours_3d * source,N_gradient_neighbours_3d * target)875 N_copy_gradient_neighbours_3d(N_gradient_neighbours_3d * source,
876 N_gradient_neighbours_3d * target)
877 {
878 int fail = 0;
879
880 G_debug(5,
881 "N_copy_gradient_neighbours_3d: copy N_gradient_neighbours_3d");
882
883 if (!source || !target)
884 return 0;
885
886 if (!(N_copy_gradient_neighbours_x(source->xt, target->xt)))
887 fail++;
888 if (!(N_copy_gradient_neighbours_x(source->xc, target->xc)))
889 fail++;
890 if (!(N_copy_gradient_neighbours_x(source->xb, target->xb)))
891 fail++;
892 if (!(N_copy_gradient_neighbours_y(source->yt, target->yt)))
893 fail++;
894 if (!(N_copy_gradient_neighbours_y(source->yc, target->yc)))
895 fail++;
896 if (!(N_copy_gradient_neighbours_y(source->yb, target->yb)))
897 fail++;
898 if (!(N_copy_gradient_neighbours_z(source->zt, target->zt)))
899 fail++;
900 if (!(N_copy_gradient_neighbours_z(source->zb, target->zb)))
901 fail++;
902
903 if (fail > 0) {
904 return 0;
905 }
906
907 return 1;
908 }
909
910 /*!
911 * \brief Allocate a N_gradient_field_2d
912 *
913 * The field arrays are of type DCELL.
914 *
915 * \param rows - number of rows of the 2d array from which the gradient should be calculated
916 * \param cols - number of cols of the 2d array from which the gradient should be calculated
917 * \return N_gradient_field_2d *
918 *
919 * */
N_alloc_gradient_field_2d(int cols,int rows)920 N_gradient_field_2d *N_alloc_gradient_field_2d(int cols, int rows)
921 {
922 N_gradient_field_2d *field;
923
924 G_debug(5,
925 "N_alloc_gradient_field_2d: allocate a N_gradient_field_2d struct");
926
927 field = (N_gradient_field_2d *) G_calloc(1, sizeof(N_gradient_field_2d));
928
929 field->x_array = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
930 field->y_array = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
931
932 field->cols = cols;
933 field->rows = rows;
934
935 return field;
936 }
937
938 /*!
939 * \brief Free's a N_gradient_neighbours_2d structure
940 *
941 * \return void
942 *
943 * */
N_free_gradient_field_2d(N_gradient_field_2d * field)944 void N_free_gradient_field_2d(N_gradient_field_2d * field)
945 {
946
947 N_free_array_2d(field->x_array);
948 N_free_array_2d(field->y_array);
949
950 G_free(field);
951
952 field = NULL;
953
954 return;
955 }
956
957 /*!
958 * \brief Copy N_gradient_field_2d structure from source to target
959 *
960 * \param source - the source N_gradient_field_2d struct
961 * \param target - the target N_gradient_field_2d struct
962 * \return int - 1 success, 0 failure while copying
963 *
964 * */
965 int
N_copy_gradient_field_2d(N_gradient_field_2d * source,N_gradient_field_2d * target)966 N_copy_gradient_field_2d(N_gradient_field_2d * source,
967 N_gradient_field_2d * target)
968 {
969 G_debug(3, "N_copy_gradient_field_2d: copy N_gradient_field_2d");
970
971 if (!source || !target)
972 return 0;
973
974 N_copy_array_2d(source->x_array, target->x_array);
975 N_copy_array_2d(source->y_array, target->y_array);
976
977 return 1;
978 }
979
980 /*! \brief Print gradient field information to stdout
981 *
982 * \param field N_gradient_2d_field *
983 * \return void
984 *
985 * */
N_print_gradient_field_2d_info(N_gradient_field_2d * field)986 void N_print_gradient_field_2d_info(N_gradient_field_2d * field)
987 {
988 fprintf(stdout, "N_gradient_field_2d \n");
989 fprintf(stdout, "Cols %i\n", field->cols);
990 fprintf(stdout, "Rows: %i\n", field->rows);
991 fprintf(stdout, "X array pointer: %p\n", field->x_array);
992 fprintf(stdout, "Y array pointer: %p\n", field->y_array);
993 fprintf(stdout, "Min %g\n", field->min);
994 fprintf(stdout, "Max %g\n", field->max);
995 fprintf(stdout, "Sum %g\n", field->sum);
996 fprintf(stdout, "Mean %g\n", field->mean);
997 fprintf(stdout, "Nonull %i\n", field->nonull);
998 fprintf(stdout, "X array info \n");
999 N_print_array_2d_info(field->x_array);
1000 fprintf(stdout, "Y array info \n");
1001 N_print_array_2d_info(field->y_array);
1002
1003 return;
1004 }
1005
1006
1007 /*!
1008 * \brief Allocate a N_gradient_field_3d
1009 *
1010 * The field arrays are always of type DCELL_TYPE.
1011 *
1012 * \param cols - number of cols of the 3d array from which the gradient should be calculated
1013 * \param rows - number of rows of the 3d array from which the gradient should be calculated
1014 * \param depths - number of depths of the 3d array from which the gradient should be calculated
1015 * \return N_gradient_field_3d *
1016 *
1017 * */
N_alloc_gradient_field_3d(int cols,int rows,int depths)1018 N_gradient_field_3d *N_alloc_gradient_field_3d(int cols, int rows, int depths)
1019 {
1020 N_gradient_field_3d *field;
1021
1022 G_debug(5,
1023 "N_alloc_gradient_field_3d: allocate a N_gradient_field_3d struct");
1024
1025 field = (N_gradient_field_3d *) G_calloc(1, sizeof(N_gradient_field_3d));
1026
1027 field->x_array = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
1028 field->y_array = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
1029 field->z_array = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
1030
1031 field->cols = cols;
1032 field->rows = rows;
1033 field->depths = depths;
1034
1035 return field;
1036 }
1037
1038
1039 /*!
1040 * \brief Free's a N_gradient_neighbours_3d structure
1041 *
1042 * \return void
1043 *
1044 * */
N_free_gradient_field_3d(N_gradient_field_3d * field)1045 void N_free_gradient_field_3d(N_gradient_field_3d * field)
1046 {
1047
1048 N_free_array_3d(field->x_array);
1049 N_free_array_3d(field->y_array);
1050 N_free_array_3d(field->z_array);
1051
1052 G_free(field);
1053
1054 field = NULL;
1055
1056 return;
1057 }
1058
1059
1060 /*!
1061 * \brief Copy N_gradient_field_3d structure from source to target
1062 *
1063 * \param source - the source N_gradient_field_3d struct
1064 * \param target - the target N_gradient_field_3d struct
1065 * \return int - 1 success, 0 failure while copying
1066 *
1067 * */
1068 int
N_copy_gradient_field_3d(N_gradient_field_3d * source,N_gradient_field_3d * target)1069 N_copy_gradient_field_3d(N_gradient_field_3d * source,
1070 N_gradient_field_3d * target)
1071 {
1072 G_debug(3, "N_copy_gradient_field_3d: copy N_gradient_field_3d");
1073
1074 if (!source || !target)
1075 return 0;
1076
1077 N_copy_array_3d(source->x_array, target->x_array);
1078 N_copy_array_3d(source->y_array, target->y_array);
1079 N_copy_array_3d(source->z_array, target->z_array);
1080
1081 return 1;
1082 }
1083
1084 /*! \brief Print gradient field information to stdout
1085 *
1086 * \param field N_gradient_3d_field *
1087 * \return void
1088 *
1089 * */
N_print_gradient_field_3d_info(N_gradient_field_3d * field)1090 void N_print_gradient_field_3d_info(N_gradient_field_3d * field)
1091 {
1092
1093 fprintf(stdout, "N_gradient_field_3d \n");
1094 fprintf(stdout, "Cols %i\n", field->cols);
1095 fprintf(stdout, "Rows: %i\n", field->rows);
1096 fprintf(stdout, "Depths %i\n", field->depths);
1097 fprintf(stdout, "X array pointer: %p\n", field->x_array);
1098 fprintf(stdout, "Y array pointer: %p\n", field->y_array);
1099 fprintf(stdout, "Z array pointer: %p\n", field->z_array);
1100 fprintf(stdout, "Min %g\n", field->min);
1101 fprintf(stdout, "Max %g\n", field->max);
1102 fprintf(stdout, "Sum %g\n", field->sum);
1103 fprintf(stdout, "Mean %g\n", field->mean);
1104 fprintf(stdout, "Nonull %i\n", field->nonull);
1105 fprintf(stdout, "X array info \n");
1106 N_print_array_3d_info(field->x_array);
1107 fprintf(stdout, "Y array info \n");
1108 N_print_array_3d_info(field->y_array);
1109 fprintf(stdout, "Z array info \n");
1110 N_print_array_3d_info(field->z_array);
1111
1112 return;
1113 }
1114