GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-d6dec75dd4
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 
34 void Rast3d_gradient_double(RASTER3D_Array_double *array, double *step,
35  RASTER3D_Array_double *grad_x,
36  RASTER3D_Array_double *grad_y,
37  RASTER3D_Array_double *grad_z)
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);
48  if (Rast_is_d_null_value(&val0))
49  Rast_set_null_value(
50  &RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth), 1,
51  DCELL_TYPE);
52  else if (Rast_is_d_null_value(&val1) || Rast_is_d_null_value(&val2))
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);
62  if (Rast_is_d_null_value(&val2))
63  Rast_set_null_value(
64  &RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth),
65  1, DCELL_TYPE);
66  else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val1))
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);
77  if (Rast_is_d_null_value(&val1))
78  Rast_set_null_value(
79  &RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth), 1,
80  DCELL_TYPE);
81  else if (Rast_is_d_null_value(&val0) ||
82  Rast_is_d_null_value(&val2))
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);
96  if (Rast_is_d_null_value(&val0))
97  Rast_set_null_value(
98  &RASTER3D_ARRAY_ACCESS(grad_y, col, 0, depth), 1,
99  DCELL_TYPE);
100  else if (Rast_is_d_null_value(&val1) || Rast_is_d_null_value(&val2))
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);
110  if (Rast_is_d_null_value(&val2))
111  Rast_set_null_value(
112  &RASTER3D_ARRAY_ACCESS(grad_y, col, array->sy - 1, depth),
113  1, DCELL_TYPE);
114  else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val1))
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);
125  if (Rast_is_d_null_value(&val1))
126  Rast_set_null_value(
127  &RASTER3D_ARRAY_ACCESS(grad_y, col, row, depth), 1,
128  DCELL_TYPE);
129  else if (Rast_is_d_null_value(&val0) ||
130  Rast_is_d_null_value(&val2))
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);
144  if (Rast_is_d_null_value(&val0))
145  Rast_set_null_value(&RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0),
146  1, DCELL_TYPE);
147  else if (Rast_is_d_null_value(&val1) || Rast_is_d_null_value(&val2))
148  RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0) = 0;
149  else
150  RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0) =
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);
157  if (Rast_is_d_null_value(&val2))
158  Rast_set_null_value(
159  &RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1), 1,
160  DCELL_TYPE);
161  else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val1))
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);
171  if (Rast_is_d_null_value(&val1))
172  Rast_set_null_value(
173  &RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth), 1,
174  DCELL_TYPE);
175  else if (Rast_is_d_null_value(&val0) ||
176  Rast_is_d_null_value(&val2))
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 }
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:98
Rast_is_d_null_value
#define Rast_is_d_null_value(dcellVal)
Definition: defs/raster.h:408
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:34
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