GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
gradient.c
Go to the documentation of this file.
1 
2 /*!
3  \file gradient.c
4 
5  \brief Gradient computation
6 
7  (C) 2014 by the GRASS Development Team
8 
9  This program is free software under the GNU General Public
10  License (>=v2). Read the file COPYING that comes with GRASS
11  for details.
12 
13  \author Anna Petrasova
14  */
15 
16 /*!
17  \brief Gradient computation
18 
19  Gradient computation (second order approximation)
20  using central differencing scheme (plus forward and backward
21  difference of second order approximation). When one or more of the cells,
22  from which the gradient for a particular cell is computed, is null,
23  gradient for that particular cell is set to 0.
24 
25  \param array pointer to RASTER3D_Array with input values
26  \param step array of x, y, z steps for gradient (resolution values)
27  \param[out] grad_x pointer to RASTER3D_Array_double with gradient in x direction
28  \param[out] grad_y pointer to RASTER3D_Array_double with gradient in y direction
29  \param[out] grad_z pointer to RASTER3D_Array_double with 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(&RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth),
50  1, DCELL_TYPE);
51  else if (Rast_is_d_null_value(&val1) || Rast_is_d_null_value(&val2))
52  RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth) = 0;
53  else
54  RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth) =
55  (-3 * val0 + 4 * val1 - val2) / (2 * step[0]);
56 
57  /* row end */
58  val0 = RASTER3D_ARRAY_ACCESS(array, array->sx - 3, row, depth);
59  val1 = RASTER3D_ARRAY_ACCESS(array, array->sx - 2, row, depth);
60  val2 = RASTER3D_ARRAY_ACCESS(array, array->sx - 1, row, depth);
61  if (Rast_is_d_null_value(&val2))
62  Rast_set_null_value(
63  &RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth),
64  1, DCELL_TYPE);
65  else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val1))
66  RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth) = 0;
67  else
68  RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth) =
69  (3 * val2 - 4 * val1 + val0) / (2 * step[0]);
70 
71  /* row */
72  for (col = 1; col < array->sx - 1; col++) {
73  val0 = RASTER3D_ARRAY_ACCESS(array, col - 1, row, depth);
74  val1 = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
75  val2 = RASTER3D_ARRAY_ACCESS(array, col + 1, row, depth);
76  if (Rast_is_d_null_value(&val1))
77  Rast_set_null_value(
78  &RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth),
79  1, DCELL_TYPE);
80  else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val2))
81  RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth) = 0;
82  else
83  RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth) =
84  (val2 - val0) / (2 * step[0]);
85  }
86  }
87  }
88  for (depth = 0; depth < array->sz; depth++) {
89  for (col = 0; col < array->sx; col++) {
90  /* col start */
91  val0 = RASTER3D_ARRAY_ACCESS(array, col, 0, depth);
92  val1 = RASTER3D_ARRAY_ACCESS(array, col, 1, depth);
93  val2 = RASTER3D_ARRAY_ACCESS(array, col, 2, depth);
94  if (Rast_is_d_null_value(&val0))
95  Rast_set_null_value(&RASTER3D_ARRAY_ACCESS(grad_y, col, 0, depth),
96  1, 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),
125  1, DCELL_TYPE);
126  else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val2))
127  RASTER3D_ARRAY_ACCESS(grad_y, col, row, depth) = 0;
128  else
129  RASTER3D_ARRAY_ACCESS(grad_y, col, row, depth) =
130  -(val2 - val0) / (2 * step[1]);
131  }
132  }
133  }
134  for (row = 0; row < array->sy; row++) {
135  for (col = 0; col < array->sx; col++) {
136  /* vertical col start */
137  val0 = RASTER3D_ARRAY_ACCESS(array, col, row, 0);
138  val1 = RASTER3D_ARRAY_ACCESS(array, col, row, 1);
139  val2 = RASTER3D_ARRAY_ACCESS(array, col, row, 2);
140  if (Rast_is_d_null_value(&val0))
141  Rast_set_null_value(&RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0),
142  1, DCELL_TYPE);
143  else if (Rast_is_d_null_value(&val1) || Rast_is_d_null_value(&val2))
144  RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0) = 0;
145  else
146  RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0) =
147  (-3 * val0 + 4 * val1 - val2) / (2 * step[2]);
148 
149  /* vertical col end */
150  val0 = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 3);
151  val1 = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 2);
152  val2 = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 1);
153  if (Rast_is_d_null_value(&val2))
154  Rast_set_null_value(
155  &RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1),
156  1, DCELL_TYPE);
157  else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val1))
158  RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1) = 0;
159  else
160  RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1) =
161  (3 * val2 - 4 * val1 + val0) / (2 * step[2]);
162  /* vertical col */
163  for (depth = 1; depth < array->sz - 1; depth++) {
164  val0 = RASTER3D_ARRAY_ACCESS(array, col, row, depth - 1);
165  val1 = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
166  val2 = RASTER3D_ARRAY_ACCESS(array, col, row, depth + 1);
167  if (Rast_is_d_null_value(&val1))
168  Rast_set_null_value(
169  &RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth),
170  1, DCELL_TYPE);
171  else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val2))
172  RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth) = 0;
173  else
174  RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth) =
175  (val2 - val0) / (2 * step[2]);
176  }
177  }
178  }
179 }
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
RASTER3D_Array_double
Definition: raster3d.h:259
Rast_is_d_null_value
#define Rast_is_d_null_value(dcellVal)
Definition: defs/raster.h:414
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_double::sz
int sz
Definition: raster3d.h:264
RASTER3D_Array_double::sx
int sx
Definition: raster3d.h:262
DCELL_TYPE
#define DCELL_TYPE
Definition: raster.h:13
RASTER3D_Array_double::sy
int sy
Definition: raster3d.h:263
RASTER3D_ARRAY_ACCESS
#define RASTER3D_ARRAY_ACCESS(arr, x, y, z)
Definition: raster3d.h:268