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 
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))
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))
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))
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))
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))
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))
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 }
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:414
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 DCELL_TYPE
Definition: raster.h:13
#define RASTER3D_ARRAY_ACCESS(arr, x, y, z)
Definition: raster3d.h:268