GRASS 8 Programmer's Manual 8.6.0dev(2026)-f6f2c534ea
Loading...
Searching...
No Matches
gradient.c
Go to the documentation of this file.
1/*!
2 \file gradient.c
3
4 \brief Gradient computation
5
6 (C) 2014 by the GRASS Development Team
7
8 This program is free software under the GNU General Public
9 License (>=v2). Read the file COPYING that comes with GRASS
10 for details.
11
12 \author Anna Petrasova
13 */
14
15/*!
16 \brief Gradient computation
17
18 Gradient computation (second order approximation)
19 using central differencing scheme (plus forward and backward
20 difference of second order approximation). When one or more of the cells,
21 from which the gradient for a particular cell is computed, is null,
22 gradient for that particular cell is set to 0.
23
24 \param array pointer to RASTER3D_Array with input values
25 \param step array of x, y, z steps for gradient (resolution values)
26 \param[out] grad_x pointer to RASTER3D_Array_double with gradient in x
27 direction \param[out] grad_y pointer to RASTER3D_Array_double with gradient
28 in y direction \param[out] grad_z pointer to RASTER3D_Array_double with
29 gradient in z direction
30
31 */
32#include <grass/raster3d.h>
33
38{
39 int col, row, depth;
40 double val0, val1, val2;
41
42 for (depth = 0; depth < array->sz; depth++) {
43 for (row = 0; row < array->sy; row++) {
44 /* row start */
45 val0 = RASTER3D_ARRAY_ACCESS(array, 0, row, depth);
46 val1 = RASTER3D_ARRAY_ACCESS(array, 1, row, depth);
47 val2 = RASTER3D_ARRAY_ACCESS(array, 2, row, depth);
50 &RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth), 1,
53 RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth) = 0;
54 else
55 RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth) =
56 (-3 * val0 + 4 * val1 - val2) / (2 * step[0]);
57
58 /* row end */
59 val0 = RASTER3D_ARRAY_ACCESS(array, array->sx - 3, row, depth);
60 val1 = RASTER3D_ARRAY_ACCESS(array, array->sx - 2, row, depth);
61 val2 = RASTER3D_ARRAY_ACCESS(array, array->sx - 1, row, depth);
64 &RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth),
65 1, DCELL_TYPE);
67 RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth) = 0;
68 else
69 RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth) =
70 (3 * val2 - 4 * val1 + val0) / (2 * step[0]);
71
72 /* row */
73 for (col = 1; col < array->sx - 1; col++) {
74 val0 = RASTER3D_ARRAY_ACCESS(array, col - 1, row, depth);
75 val1 = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
76 val2 = RASTER3D_ARRAY_ACCESS(array, col + 1, row, depth);
79 &RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth), 1,
81 else if (Rast_is_d_null_value(&val0) ||
83 RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth) = 0;
84 else
85 RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth) =
86 (val2 - val0) / (2 * step[0]);
87 }
88 }
89 }
90 for (depth = 0; depth < array->sz; depth++) {
91 for (col = 0; col < array->sx; col++) {
92 /* col start */
93 val0 = RASTER3D_ARRAY_ACCESS(array, col, 0, depth);
94 val1 = RASTER3D_ARRAY_ACCESS(array, col, 1, depth);
95 val2 = RASTER3D_ARRAY_ACCESS(array, col, 2, depth);
98 &RASTER3D_ARRAY_ACCESS(grad_y, col, 0, depth), 1,
101 RASTER3D_ARRAY_ACCESS(grad_y, col, 0, depth) = 0;
102 else
103 RASTER3D_ARRAY_ACCESS(grad_y, col, 0, depth) =
104 -(-3 * val0 + 4 * val1 - val2) / (2 * step[1]);
105
106 /* col end */
107 val0 = RASTER3D_ARRAY_ACCESS(array, col, array->sy - 3, depth);
108 val1 = RASTER3D_ARRAY_ACCESS(array, col, array->sy - 2, depth);
109 val2 = RASTER3D_ARRAY_ACCESS(array, col, array->sy - 1, depth);
112 &RASTER3D_ARRAY_ACCESS(grad_y, col, array->sy - 1, depth),
113 1, DCELL_TYPE);
115 RASTER3D_ARRAY_ACCESS(grad_y, col, array->sy - 1, depth) = 0;
116 else
117 RASTER3D_ARRAY_ACCESS(grad_y, col, array->sy - 1, depth) =
118 -(3 * val2 - 4 * val1 + val0) / (2 * step[1]);
119
120 /* col */
121 for (row = 1; row < array->sy - 1; row++) {
122 val0 = RASTER3D_ARRAY_ACCESS(array, col, row - 1, depth);
123 val1 = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
124 val2 = RASTER3D_ARRAY_ACCESS(array, col, row + 1, depth);
127 &RASTER3D_ARRAY_ACCESS(grad_y, col, row, depth), 1,
128 DCELL_TYPE);
129 else if (Rast_is_d_null_value(&val0) ||
131 RASTER3D_ARRAY_ACCESS(grad_y, col, row, depth) = 0;
132 else
133 RASTER3D_ARRAY_ACCESS(grad_y, col, row, depth) =
134 -(val2 - val0) / (2 * step[1]);
135 }
136 }
137 }
138 for (row = 0; row < array->sy; row++) {
139 for (col = 0; col < array->sx; col++) {
140 /* vertical col start */
141 val0 = RASTER3D_ARRAY_ACCESS(array, col, row, 0);
142 val1 = RASTER3D_ARRAY_ACCESS(array, col, row, 1);
143 val2 = RASTER3D_ARRAY_ACCESS(array, col, row, 2);
146 1, DCELL_TYPE);
148 RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0) = 0;
149 else
151 (-3 * val0 + 4 * val1 - val2) / (2 * step[2]);
152
153 /* vertical col end */
154 val0 = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 3);
155 val1 = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 2);
156 val2 = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 1);
159 &RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1), 1,
160 DCELL_TYPE);
162 RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1) = 0;
163 else
164 RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1) =
165 (3 * val2 - 4 * val1 + val0) / (2 * step[2]);
166 /* vertical col */
167 for (depth = 1; depth < array->sz - 1; depth++) {
168 val0 = RASTER3D_ARRAY_ACCESS(array, col, row, depth - 1);
169 val1 = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
170 val2 = RASTER3D_ARRAY_ACCESS(array, col, row, depth + 1);
173 &RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth), 1,
174 DCELL_TYPE);
175 else if (Rast_is_d_null_value(&val0) ||
177 RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth) = 0;
178 else
179 RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth) =
180 (val2 - val0) / (2 * step[2]);
181 }
182 }
183 }
184}
void Rast_set_null_value(void *, int, RASTER_MAP_TYPE)
To set one or more raster values to null.
Definition null_val.c:98
#define Rast_is_d_null_value(dcellVal)
void Rast3d_gradient_double(RASTER3D_Array_double *array, double *step, RASTER3D_Array_double *grad_x, RASTER3D_Array_double *grad_y, RASTER3D_Array_double *grad_z)
Gradient computation.
Definition gradient.c:34
#define RASTER3D_ARRAY_ACCESS(arr, x, y, z)
Definition raster3d.h:267
#define DCELL_TYPE
Definition raster.h:13