GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-87d2d42d3b
Go to the documentation of this file.
1 /*!
3
5
6  (C) 2014 by the GRASS Development Team
7
8  This program is free software under the GNU General Public
10  for details.
11
12  \author Anna Petrasova
13  */
14
15 /*!
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)
28  in y direction \param[out] grad_z pointer to RASTER3D_Array_double with
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);
48  if (Rast_is_d_null_value(&val0))
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))
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))
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))
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))
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))
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))
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))
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))
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 }
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)
Definition: defs/raster.h:405