GRASS 8 Programmer's Manual 8.6.0dev(2026)-4bb960b182
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 SPDX-FileCopyrightText: 2014 GRASS Development Team
7 SPDX-License-Identifier: GPL-2.0-or-later
8
9 \author Anna Petrasova
10 */
11
12/*!
13 \brief Gradient computation
14
15 Gradient computation (second order approximation)
16 using central differencing scheme (plus forward and backward
17 difference of second order approximation). When one or more of the cells,
18 from which the gradient for a particular cell is computed, is null,
19 gradient for that particular cell is set to 0.
20
21 \param array pointer to RASTER3D_Array with input values
22 \param step array of x, y, z steps for gradient (resolution values)
23 \param[out] grad_x pointer to RASTER3D_Array_double with gradient in x
24 direction \param[out] grad_y pointer to RASTER3D_Array_double with gradient
25 in y direction \param[out] grad_z pointer to RASTER3D_Array_double with
26 gradient in z direction
27
28 */
29#include <grass/raster3d.h>
30
31void Rast3d_gradient_double(RASTER3D_Array_double *array, double *step,
32 RASTER3D_Array_double *grad_x,
33 RASTER3D_Array_double *grad_y,
34 RASTER3D_Array_double *grad_z)
35{
36 int col, row, depth;
37 double val0, val1, val2;
38
39 for (depth = 0; depth < array->sz; depth++) {
40 for (row = 0; row < array->sy; row++) {
41 /* row start */
42 val0 = RASTER3D_ARRAY_ACCESS(array, 0, row, depth);
43 val1 = RASTER3D_ARRAY_ACCESS(array, 1, row, depth);
44 val2 = RASTER3D_ARRAY_ACCESS(array, 2, row, depth);
45 if (Rast_is_d_null_value(&val0))
46 Rast_set_null_value(
47 &RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth), 1,
48 DCELL_TYPE);
49 else if (Rast_is_d_null_value(&val1) || Rast_is_d_null_value(&val2))
50 RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth) = 0;
51 else
52 RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth) =
53 (-3 * val0 + 4 * val1 - val2) / (2 * step[0]);
54
55 /* row end */
56 val0 = RASTER3D_ARRAY_ACCESS(array, array->sx - 3, row, depth);
57 val1 = RASTER3D_ARRAY_ACCESS(array, array->sx - 2, row, depth);
58 val2 = RASTER3D_ARRAY_ACCESS(array, array->sx - 1, row, depth);
59 if (Rast_is_d_null_value(&val2))
60 Rast_set_null_value(
61 &RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth),
62 1, DCELL_TYPE);
63 else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val1))
64 RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth) = 0;
65 else
66 RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth) =
67 (3 * val2 - 4 * val1 + val0) / (2 * step[0]);
68
69 /* row */
70 for (col = 1; col < array->sx - 1; col++) {
71 val0 = RASTER3D_ARRAY_ACCESS(array, col - 1, row, depth);
72 val1 = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
73 val2 = RASTER3D_ARRAY_ACCESS(array, col + 1, row, depth);
74 if (Rast_is_d_null_value(&val1))
75 Rast_set_null_value(
76 &RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth), 1,
77 DCELL_TYPE);
78 else if (Rast_is_d_null_value(&val0) ||
79 Rast_is_d_null_value(&val2))
80 RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth) = 0;
81 else
82 RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth) =
83 (val2 - val0) / (2 * step[0]);
84 }
85 }
86 }
87 for (depth = 0; depth < array->sz; depth++) {
88 for (col = 0; col < array->sx; col++) {
89 /* col start */
90 val0 = RASTER3D_ARRAY_ACCESS(array, col, 0, depth);
91 val1 = RASTER3D_ARRAY_ACCESS(array, col, 1, depth);
92 val2 = RASTER3D_ARRAY_ACCESS(array, col, 2, depth);
93 if (Rast_is_d_null_value(&val0))
94 Rast_set_null_value(
95 &RASTER3D_ARRAY_ACCESS(grad_y, col, 0, depth), 1,
96 DCELL_TYPE);
97 else if (Rast_is_d_null_value(&val1) || Rast_is_d_null_value(&val2))
98 RASTER3D_ARRAY_ACCESS(grad_y, col, 0, depth) = 0;
99 else
100 RASTER3D_ARRAY_ACCESS(grad_y, col, 0, depth) =
101 -(-3 * val0 + 4 * val1 - val2) / (2 * step[1]);
102
103 /* col end */
104 val0 = RASTER3D_ARRAY_ACCESS(array, col, array->sy - 3, depth);
105 val1 = RASTER3D_ARRAY_ACCESS(array, col, array->sy - 2, depth);
106 val2 = RASTER3D_ARRAY_ACCESS(array, col, array->sy - 1, depth);
107 if (Rast_is_d_null_value(&val2))
108 Rast_set_null_value(
109 &RASTER3D_ARRAY_ACCESS(grad_y, col, array->sy - 1, depth),
110 1, DCELL_TYPE);
111 else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val1))
112 RASTER3D_ARRAY_ACCESS(grad_y, col, array->sy - 1, depth) = 0;
113 else
114 RASTER3D_ARRAY_ACCESS(grad_y, col, array->sy - 1, depth) =
115 -(3 * val2 - 4 * val1 + val0) / (2 * step[1]);
116
117 /* col */
118 for (row = 1; row < array->sy - 1; row++) {
119 val0 = RASTER3D_ARRAY_ACCESS(array, col, row - 1, depth);
120 val1 = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
121 val2 = RASTER3D_ARRAY_ACCESS(array, col, row + 1, depth);
122 if (Rast_is_d_null_value(&val1))
123 Rast_set_null_value(
124 &RASTER3D_ARRAY_ACCESS(grad_y, col, row, depth), 1,
125 DCELL_TYPE);
126 else if (Rast_is_d_null_value(&val0) ||
127 Rast_is_d_null_value(&val2))
128 RASTER3D_ARRAY_ACCESS(grad_y, col, row, depth) = 0;
129 else
130 RASTER3D_ARRAY_ACCESS(grad_y, col, row, depth) =
131 -(val2 - val0) / (2 * step[1]);
132 }
133 }
134 }
135 for (row = 0; row < array->sy; row++) {
136 for (col = 0; col < array->sx; col++) {
137 /* vertical col start */
138 val0 = RASTER3D_ARRAY_ACCESS(array, col, row, 0);
139 val1 = RASTER3D_ARRAY_ACCESS(array, col, row, 1);
140 val2 = RASTER3D_ARRAY_ACCESS(array, col, row, 2);
141 if (Rast_is_d_null_value(&val0))
142 Rast_set_null_value(&RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0),
143 1, DCELL_TYPE);
144 else if (Rast_is_d_null_value(&val1) || Rast_is_d_null_value(&val2))
145 RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0) = 0;
146 else
147 RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0) =
148 (-3 * val0 + 4 * val1 - val2) / (2 * step[2]);
149
150 /* vertical col end */
151 val0 = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 3);
152 val1 = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 2);
153 val2 = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 1);
154 if (Rast_is_d_null_value(&val2))
155 Rast_set_null_value(
156 &RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1), 1,
157 DCELL_TYPE);
158 else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val1))
159 RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1) = 0;
160 else
161 RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1) =
162 (3 * val2 - 4 * val1 + val0) / (2 * step[2]);
163 /* vertical col */
164 for (depth = 1; depth < array->sz - 1; depth++) {
165 val0 = RASTER3D_ARRAY_ACCESS(array, col, row, depth - 1);
166 val1 = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
167 val2 = RASTER3D_ARRAY_ACCESS(array, col, row, depth + 1);
168 if (Rast_is_d_null_value(&val1))
169 Rast_set_null_value(
170 &RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth), 1,
171 DCELL_TYPE);
172 else if (Rast_is_d_null_value(&val0) ||
173 Rast_is_d_null_value(&val2))
174 RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth) = 0;
175 else
176 RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth) =
177 (val2 - val0) / (2 * step[2]);
178 }
179 }
180 }
181}
AMI_STREAM
Definition ami_stream.h:153
Rast_set_null_value
void Rast_set_null_value(void *, int, RASTER_MAP_TYPE)
To set one or more raster values to null.
Definition null_val.c:96
Rast_is_d_null_value
#define Rast_is_d_null_value(dcellVal)
Definition defs/raster.h:417
Rast3d_gradient_double
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:31
RASTER3D_ARRAY_ACCESS
#define RASTER3D_ARRAY_ACCESS(arr, x, y, z)
Definition raster3d.h:267
DCELL_TYPE
#define DCELL_TYPE
Definition raster.h:13
RASTER3D_Array_double
Definition raster3d.h:259
RASTER3D_Array_double::sx
int sx
Definition raster3d.h:261
RASTER3D_Array_double::sz
int sz
Definition raster3d.h:263
RASTER3D_Array_double::sy
int sy
Definition raster3d.h:262