GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-535c39c9fc
sample.c
Go to the documentation of this file.
1 /*!
2  \file lib/raster/sample.c
3 
4  \brief Raster library - Sampling methods (extract a cell value from
5  raster map)
6 
7  1/2006: moved to libgis from v.sample/v.drape for clone removal
8 
9  (C) 2001-2009 by the GRASS Development Team
10 
11  This program is free software under the GNU General Public License
12  (>=v2). Read the file COPYING that comes with GRASS for details.
13 
14  \author James Darrell McCauley <darrell mccauley-usa.com>,
15  http://mccauley-usa.com/
16  */
17 
18 #include <string.h>
19 #include <unistd.h>
20 #include <math.h>
21 
22 #include <grass/gis.h>
23 #include <grass/raster.h>
24 #include <grass/glocale.h>
25 
26 /* prototypes */
27 static double scancatlabel(const char *);
28 
29 /*!
30  * \brief Extract a cell value from raster map.
31  *
32  * Extract a cell value from raster map at given northing and easting
33  * with a sampled 3x3 window using a specified interpolation method.
34  *
35  * - NEAREST neighbor interpolation
36  * - BILINEAR bilinear interpolation
37  * - CUBIC cubic interpolation
38  *
39  * \param fd file descriptor
40  * \param window region settings
41  * \param cats categories
42  * \param north northing position
43  * \param east easting position
44  * \param usedesc flag to scan category label
45  * \param itype interpolation method
46  *
47  * \return cell value at given position
48  */
49 DCELL Rast_get_sample(int fd, const struct Cell_head *window,
50  struct Categories *cats, double north, double east,
51  int usedesc, INTERP_TYPE itype)
52 {
53  double retval;
54 
55  switch (itype) {
56  case INTERP_NEAREST:
57  retval =
58  Rast_get_sample_nearest(fd, window, cats, north, east, usedesc);
59  break;
60  case INTERP_BILINEAR:
61  retval =
62  Rast_get_sample_bilinear(fd, window, cats, north, east, usedesc);
63  break;
64  case INTERP_BICUBIC:
65  retval = Rast_get_sample_cubic(fd, window, cats, north, east, usedesc);
66  break;
67  default:
68  G_fatal_error("Rast_get_sample: %s", _("Unknown interpolation type"));
69  }
70 
71  return retval;
72 }
73 
74 /*!
75  * \brief Extract a cell value from raster map (neighbor interpolation)
76  *
77  * Extract a cell value from raster map at given northing and easting
78  * with a sampled 3x3 window using a neighbor interpolation.
79  *
80  * \param fd file descriptor
81  * \param window region settings
82  * \param cats categories
83  * \param north northing position
84  * \param east easting position
85  * \param usedesc flag to scan category label
86  *
87  * \return cell value at given position
88  */
89 DCELL Rast_get_sample_nearest(int fd, const struct Cell_head *window,
90  struct Categories *cats, double north,
91  double east, int usedesc)
92 {
93  int row, col;
94  DCELL result;
95  DCELL *maprow = Rast_allocate_d_buf();
96 
97  /* convert northing and easting to row and col, resp */
98  row = (int)floor(Rast_northing_to_row(north, window));
99  col = (int)floor(Rast_easting_to_col(east, window));
100 
101  if (row < 0 || row >= Rast_window_rows() || col < 0 ||
102  col >= Rast_window_cols()) {
103  Rast_set_d_null_value(&result, 1);
104  goto done;
105  }
106 
107  Rast_get_d_row(fd, maprow, row);
108 
109  if (Rast_is_d_null_value(&maprow[col])) {
110  Rast_set_d_null_value(&result, 1);
111  goto done;
112  }
113 
114  if (usedesc) {
115  char *buf = Rast_get_c_cat((CELL *)&(maprow[col]), cats);
116 
117  G_squeeze(buf);
118  result = scancatlabel(buf);
119  }
120  else
121  result = maprow[col];
122 
123 done:
124  G_free(maprow);
125 
126  return result;
127 }
128 
129 /*!
130  * \brief Extract a cell value from raster map (bilinear interpolation).
131  *
132  * Extract a cell value from raster map at given northing and easting
133  * with a sampled 3x3 window using a bilinear interpolation.
134  *
135  * \param fd file descriptor
136  * \param window region settings
137  * \param cats categories
138  * \param north northing position
139  * \param east easting position
140  * \param usedesc flag to scan category label
141  *
142  * \return cell value at given position
143  */
144 DCELL Rast_get_sample_bilinear(int fd, const struct Cell_head *window,
145  struct Categories *cats, double north,
146  double east, int usedesc)
147 {
148  int row, col;
149  double grid[2][2];
150  DCELL *arow = Rast_allocate_d_buf();
151  DCELL *brow = Rast_allocate_d_buf();
152  double frow, fcol, trow, tcol;
153  DCELL result;
154 
155  frow = Rast_northing_to_row(north, window);
156  fcol = Rast_easting_to_col(east, window);
157 
158  /* convert northing and easting to row and col, resp */
159  row = (int)floor(frow - 0.5);
160  col = (int)floor(fcol - 0.5);
161 
162  trow = frow - row - 0.5;
163  tcol = fcol - col - 0.5;
164 
165  if (row < 0 || row + 1 >= Rast_window_rows() || col < 0 ||
166  col + 1 >= Rast_window_cols()) {
167  Rast_set_d_null_value(&result, 1);
168  goto done;
169  }
170 
171  Rast_get_d_row(fd, arow, row);
172  Rast_get_d_row(fd, brow, row + 1);
173 
174  if (Rast_is_d_null_value(&arow[col]) ||
175  Rast_is_d_null_value(&arow[col + 1]) ||
176  Rast_is_d_null_value(&brow[col]) ||
177  Rast_is_d_null_value(&brow[col + 1])) {
178  Rast_set_d_null_value(&result, 1);
179  goto done;
180  }
181 
182  /*-
183  * now were ready to do bilinear interpolation over
184  * arow[col], arow[col+1],
185  * brow[col], brow[col+1]
186  */
187 
188  if (usedesc) {
189  char *buf;
190 
191  G_squeeze(buf = Rast_get_c_cat((int *)&(arow[col]), cats));
192  grid[0][0] = scancatlabel(buf);
193  G_squeeze(buf = Rast_get_c_cat((CELL *)&(arow[col + 1]), cats));
194  grid[0][1] = scancatlabel(buf);
195  G_squeeze(buf = Rast_get_c_cat((CELL *)&(brow[col]), cats));
196  grid[1][0] = scancatlabel(buf);
197  G_squeeze(buf = Rast_get_c_cat((CELL *)&(brow[col + 1]), cats));
198  grid[1][1] = scancatlabel(buf);
199  }
200  else {
201  grid[0][0] = arow[col];
202  grid[0][1] = arow[col + 1];
203  grid[1][0] = brow[col];
204  grid[1][1] = brow[col + 1];
205  }
206 
207  result = Rast_interp_bilinear(tcol, trow, grid[0][0], grid[0][1],
208  grid[1][0], grid[1][1]);
209 
210 done:
211  G_free(arow);
212  G_free(brow);
213 
214  return result;
215 }
216 
217 /*!
218  * \brief Extract a cell value from raster map (cubic interpolation).
219  *
220  * Extract a cell value from raster map at given northing and easting
221  * with a sampled 3x3 window using a cubic interpolation.
222  *
223  * \param fd file descriptor
224  * \param window region settings
225  * \param cats categories
226  * \param north northing position
227  * \param east easting position
228  * \param usedesc flag to scan category label
229  *
230  * \return cell value at given position
231  */
232 DCELL Rast_get_sample_cubic(int fd, const struct Cell_head *window,
233  struct Categories *cats, double north, double east,
234  int usedesc)
235 {
236  int i, j, row, col;
237  double grid[4][4];
238  DCELL *rows[4];
239  double frow, fcol, trow, tcol;
240  DCELL result;
241 
242  for (i = 0; i < 4; i++)
243  rows[i] = Rast_allocate_d_buf();
244 
245  frow = Rast_northing_to_row(north, window);
246  fcol = Rast_easting_to_col(east, window);
247 
248  /* convert northing and easting to row and col, resp */
249  row = (int)floor(frow - 1.5);
250  col = (int)floor(fcol - 1.5);
251 
252  trow = frow - row - 1.5;
253  tcol = fcol - col - 1.5;
254 
255  if (row < 0 || row + 3 >= Rast_window_rows() || col < 0 ||
256  col + 3 >= Rast_window_cols()) {
257  Rast_set_d_null_value(&result, 1);
258  goto done;
259  }
260 
261  for (i = 0; i < 4; i++)
262  Rast_get_d_row(fd, rows[i], row + i);
263 
264  for (i = 0; i < 4; i++)
265  for (j = 0; j < 4; j++)
266  if (Rast_is_d_null_value(&rows[i][col + j])) {
267  Rast_set_d_null_value(&result, 1);
268  goto done;
269  }
270 
271  /*
272  * now were ready to do cubic interpolation over
273  * arow[col], arow[col+1], arow[col+2], arow[col+3],
274  * brow[col], brow[col+1], brow[col+2], brow[col+3],
275  * crow[col], crow[col+1], crow[col+2], crow[col+3],
276  * drow[col], drow[col+1], drow[col+2], drow[col+3],
277  */
278 
279  if (usedesc) {
280  char *buf;
281 
282  for (i = 0; i < 4; i++) {
283  for (j = 0; j < 4; j++) {
284  G_squeeze(
285  buf = Rast_get_c_cat((CELL *)&(rows[i][col + j]), cats));
286  grid[i][j] = scancatlabel(buf);
287  }
288  }
289  }
290  else {
291  for (i = 0; i < 4; i++)
292  for (j = 0; j < 4; j++)
293  grid[i][j] = rows[i][col + j];
294  }
295 
296  result = Rast_interp_bicubic(
297  tcol, trow, grid[0][0], grid[0][1], grid[0][2], grid[0][3], grid[1][0],
298  grid[1][1], grid[1][2], grid[1][3], grid[2][0], grid[2][1], grid[2][2],
299  grid[2][3], grid[3][0], grid[3][1], grid[3][2], grid[3][3]);
300 
301 done:
302  for (i = 0; i < 4; i++)
303  G_free(rows[i]);
304 
305  return result;
306 }
307 
308 static double scancatlabel(const char *str)
309 {
310  double val;
311 
312  if (strcmp(str, "no data") != 0)
313  sscanf(str, "%lf", &val);
314  else {
315  G_warning(_("\"no data\" label found; setting to zero"));
316  val = 0.0;
317  }
318 
319  return val;
320 }
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
void G_squeeze(char *)
Remove superfluous white space.
Definition: strings.c:446
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
DCELL Rast_interp_bilinear(double, double, DCELL, DCELL, DCELL, DCELL)
Definition: interp.c:26
double Rast_easting_to_col(double, const struct Cell_head *)
Easting to column.
DCELL Rast_interp_bicubic(double, double, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL)
Definition: interp.c:44
void Rast_set_d_null_value(DCELL *, int)
To set a number of DCELL raster values to NULL.
Definition: null_val.c:153
void Rast_get_d_row(int, DCELL *, int)
Get raster row (DCELL type)
DCELL * Rast_allocate_d_buf(void)
Allocates memory for a raster map of type DCELL.
Definition: alloc_cell.c:107
char * Rast_get_c_cat(CELL *, struct Categories *)
Get a raster category label (CELL)
Definition: raster/cats.c:321
int Rast_window_cols(void)
Number of columns in active window.
int Rast_window_rows(void)
Number of rows in active window.
Definition: raster/window.c:87
double Rast_northing_to_row(double, const struct Cell_head *)
Northing to row.
#define Rast_is_d_null_value(dcellVal)
Definition: defs/raster.h:405
double DCELL
Definition: gis.h:626
int CELL
Definition: gis.h:625
#define _(str)
Definition: glocale.h:10
#define INTERP_BILINEAR
Definition: raster.h:21
#define INTERP_NEAREST
Definition: raster.h:20
#define INTERP_BICUBIC
Definition: raster.h:22
int INTERP_TYPE
Definition: raster.h:28
DCELL Rast_get_sample_cubic(int fd, const struct Cell_head *window, struct Categories *cats, double north, double east, int usedesc)
Extract a cell value from raster map (cubic interpolation).
Definition: sample.c:232
DCELL Rast_get_sample_nearest(int fd, const struct Cell_head *window, struct Categories *cats, double north, double east, int usedesc)
Extract a cell value from raster map (neighbor interpolation)
Definition: sample.c:89
DCELL Rast_get_sample_bilinear(int fd, const struct Cell_head *window, struct Categories *cats, double north, double east, int usedesc)
Extract a cell value from raster map (bilinear interpolation).
Definition: sample.c:144
DCELL Rast_get_sample(int fd, const struct Cell_head *window, struct Categories *cats, double north, double east, int usedesc, INTERP_TYPE itype)
Extract a cell value from raster map.
Definition: sample.c:49
2D/3D raster map header (used also for region)
Definition: gis.h:437