GRASS GIS 8 Programmer's Manual  8.3.2(2024)-exported
lidar/raster.c
Go to the documentation of this file.
1 #include <grass/config.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include <grass/lidar.h>
6 
7 /*------------------------------------------------------------------------------------------------*/
8 void P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
9  struct bound_box General, struct bound_box Overlap,
10  double **obs, double *param, int *line_num, double pe,
11  double pn, double overlap, int nsplx, int nsply,
12  int num_points, int bilin, struct line_cats *categories,
13  dbDriver *driver, double mean, char *tab_name)
14 {
15  int i;
16  char buf[1024];
17  dbString sql;
18 
19  double interpolation, csi, eta, weight;
20  struct line_pnts *point;
21 
22  point = Vect_new_line_struct();
23 
24  db_begin_transaction(driver);
25 
26  for (i = 0; i < num_points; i++) {
27 
28  if (Vect_point_in_box(obs[i][0], obs[i][1], mean,
29  &General)) { /*Here mean is just for asking if obs
30  point is in box */
31 
32  if (bilin)
33  interpolation = dataInterpolateBilin(
34  obs[i][0], obs[i][1], pe, pn, nsplx, nsply,
35  Elaboration->west, Elaboration->south, param);
36  else
37  interpolation = dataInterpolateBicubic(
38  obs[i][0], obs[i][1], pe, pn, nsplx, nsply,
39  Elaboration->west, Elaboration->south, param);
40 
41  interpolation += mean;
42  Vect_copy_xyz_to_pnts(point, &obs[i][0], &obs[i][1], &interpolation,
43  1);
44 
45  if (Vect_point_in_box(obs[i][0], obs[i][1], interpolation,
46  &Overlap)) { /*(5) */
47  Vect_write_line(Out, GV_POINT, point, categories);
48  }
49  else {
50  db_init_string(&sql);
51 
52  sprintf(buf, "INSERT INTO %s (ID, X, Y, Interp)", tab_name);
53  db_append_string(&sql, buf);
54 
55  sprintf(buf, " VALUES (");
56  db_append_string(&sql, buf);
57  sprintf(buf, "%d, %f, %f, ", line_num[i], obs[i][0], obs[i][1]);
58  db_append_string(&sql, buf);
59 
60  if ((*point->x > Overlap.E) && (*point->x < General.E)) {
61  if ((*point->y > Overlap.N) &&
62  (*point->y < General.N)) { /*(3) */
63  csi = (General.E - *point->x) / overlap;
64  eta = (General.N - *point->y) / overlap;
65  weight = csi * eta;
66  *point->z = weight * interpolation;
67 
68  sprintf(buf, "%lf", *point->z);
69  db_append_string(&sql, buf);
70  sprintf(buf, ")");
71  db_append_string(&sql, buf);
72 
73  if (db_execute_immediate(driver, &sql) != DB_OK)
74  G_fatal_error(_("Unable to access table <%s>"),
75  tab_name);
76  }
77  else if ((*point->y < Overlap.S) &&
78  (*point->y > General.S)) { /*(1) */
79  csi = (General.E - *point->x) / overlap;
80  eta = (*point->y - General.S) / overlap;
81  weight = csi * eta;
82  *point->z = weight * interpolation;
83 
84  sprintf(buf, "%lf", *point->z);
85  db_append_string(&sql, buf);
86  sprintf(buf, ")");
87  db_append_string(&sql, buf);
88 
89  if (db_execute_immediate(driver, &sql) != DB_OK)
90  G_fatal_error(_("Unable to access table <%s>"),
91  tab_name);
92  }
93  else if ((*point->y <= Overlap.N) &&
94  (*point->y >= Overlap.S)) { /*(1) */
95  weight = (General.E - *point->x) / overlap;
96  *point->z = weight * interpolation;
97 
98  sprintf(buf, "%lf", *point->z);
99  db_append_string(&sql, buf);
100  sprintf(buf, ")");
101  db_append_string(&sql, buf);
102 
103  if (db_execute_immediate(driver, &sql) != DB_OK)
104  G_fatal_error(_("Unable to access table <%s>"),
105  tab_name);
106  }
107  }
108  else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
109  if ((*point->y > Overlap.N) &&
110  (*point->y < General.N)) { /*(4) */
111  csi = (*point->x - General.W) / overlap;
112  eta = (General.N - *point->y) / overlap;
113  weight = eta * csi;
114  *point->z = weight * interpolation;
115 
116  sprintf(buf, "%lf", *point->z);
117  db_append_string(&sql, buf);
118  sprintf(buf, ")");
119  db_append_string(&sql, buf);
120 
121  if (db_execute_immediate(driver, &sql) != DB_OK)
122  G_fatal_error(_("Unable to access table <%s>"),
123  tab_name);
124  }
125  else if ((*point->y < Overlap.S) &&
126  (*point->y > General.S)) { /*(2) */
127  csi = (*point->x - General.W) / overlap;
128  eta = (*point->y - General.S) / overlap;
129  weight = csi * eta;
130  *point->z = weight * interpolation;
131 
132  sprintf(buf, "%lf", *point->z);
133  db_append_string(&sql, buf);
134  sprintf(buf, ")");
135  db_append_string(&sql, buf);
136 
137  if (db_execute_immediate(driver, &sql) != DB_OK)
138  G_fatal_error(_("Unable to access table <%s>"),
139  tab_name);
140  }
141  else if ((*point->y >= Overlap.S) &&
142  (*point->y <= Overlap.N)) { /*(2) */
143  weight = (*point->x - General.W) / overlap;
144  *point->z = weight * interpolation;
145 
146  sprintf(buf, "%lf", *point->z);
147  db_append_string(&sql, buf);
148  sprintf(buf, ")");
149  db_append_string(&sql, buf);
150 
151  if (db_execute_immediate(driver, &sql) != DB_OK)
152  G_fatal_error(_("Unable to access table <%s>"),
153  tab_name);
154  }
155  }
156  else if ((*point->x >= Overlap.W) && (*point->x <= Overlap.E)) {
157  if ((*point->y > Overlap.N) &&
158  (*point->y < General.N)) { /*(3) */
159  weight = (General.N - *point->y) / overlap;
160  *point->z = weight * interpolation;
161 
162  sprintf(buf, "%lf", *point->z);
163  db_append_string(&sql, buf);
164  sprintf(buf, ")");
165  db_append_string(&sql, buf);
166 
167  if (db_execute_immediate(driver, &sql) != DB_OK)
168  G_fatal_error(_("Unable to access table <%s>"),
169  tab_name);
170  }
171  else if ((*point->y < Overlap.S) &&
172  (*point->y > General.S)) { /*(1) */
173  weight = (*point->y - General.S) / overlap;
174  *point->z = (1 - weight) * interpolation;
175 
176  sprintf(buf, "%lf", *point->z);
177  db_append_string(&sql, buf);
178  sprintf(buf, ")");
179  db_append_string(&sql, buf);
180 
181  if (db_execute_immediate(driver, &sql) != DB_OK)
182  G_fatal_error(_("Unable to access table <%s>"),
183  tab_name);
184  }
185  }
186  }
187  }
188  /*IF*/}
189  /*FOR*/ db_commit_transaction(driver);
190 
191  return;
192 }
193 
194 /*------------------------------------------------------------------------------------------------*/
195 int P_Regular_Points(struct Cell_head *Elaboration, struct Cell_head *Original,
196  struct bound_box General, struct bound_box Overlap,
197  SEGMENT *out_seg, double *param, double passoN,
198  double passoE, double overlap, double mean, int nsplx,
199  int nsply, int nrows, int ncols, int bilin)
200 {
201 
202  int col, row, startcol, endcol, startrow, endrow;
203  double X, Y, interpolation, weight, csi, eta, dval;
204 
205  /* G_get_window(&Original); */
206  if (Original->north > General.N)
207  startrow = (Original->north - General.N) / Original->ns_res - 1;
208  else
209  startrow = 0;
210  if (Original->north > General.S) {
211  endrow = (Original->north - General.S) / Original->ns_res + 1;
212  if (endrow > nrows)
213  endrow = nrows;
214  }
215  else
216  endrow = nrows;
217  if (General.W > Original->west)
218  startcol = (General.W - Original->west) / Original->ew_res - 1;
219  else
220  startcol = 0;
221  if (General.E > Original->west) {
222  endcol = (General.E - Original->west) / Original->ew_res + 1;
223  if (endcol > ncols)
224  endcol = ncols;
225  }
226  else
227  endcol = ncols;
228 
229  for (row = startrow; row < endrow; row++) {
230  for (col = startcol; col < endcol; col++) {
231 
232  X = Rast_col_to_easting((double)(col) + 0.5, Original);
233  Y = Rast_row_to_northing((double)(row) + 0.5, Original);
234 
235  if (Vect_point_in_box(X, Y, mean,
236  &General)) { /* Here, mean is just for asking
237  if obs point is in box */
238 
239  if (bilin)
240  interpolation = dataInterpolateBilin(
241  X, Y, passoE, passoN, nsplx, nsply, Elaboration->west,
242  Elaboration->south, param);
243  else
244  interpolation = dataInterpolateBicubic(
245  X, Y, passoE, passoN, nsplx, nsply, Elaboration->west,
246  Elaboration->south, param);
247 
248  interpolation += mean;
249 
250  if (Vect_point_in_box(X, Y, interpolation,
251  &Overlap)) { /* (5) */
252  dval = interpolation;
253  }
254  else {
255  Segment_get(out_seg, &dval, row, col);
256  if ((X > Overlap.E) && (X < General.E)) {
257  if ((Y > Overlap.N) && (Y < General.N)) { /* (3) */
258  csi = (General.E - X) / overlap;
259  eta = (General.N - Y) / overlap;
260  weight = csi * eta;
261  interpolation *= weight;
262  dval += interpolation;
263  }
264  else if ((Y < Overlap.S) && (Y > General.S)) { /* (1) */
265  csi = (General.E - X) / overlap;
266  eta = (Y - General.S) / overlap;
267  weight = csi * eta;
268  interpolation *= weight;
269  dval = interpolation;
270  }
271  else if ((Y >= Overlap.S) &&
272  (Y <= Overlap.N)) { /* (1) */
273  weight = (General.E - X) / overlap;
274  interpolation *= weight;
275  dval = interpolation;
276  }
277  }
278  else if ((X < Overlap.W) && (X > General.W)) {
279  if ((Y > Overlap.N) && (Y < General.N)) { /* (4) */
280  csi = (X - General.W) / overlap;
281  eta = (General.N - Y) / overlap;
282  weight = eta * csi;
283  interpolation *= weight;
284  dval += interpolation;
285  }
286  else if ((Y < Overlap.S) && (Y > General.S)) { /* (2) */
287  csi = (X - General.W) / overlap;
288  eta = (Y - General.S) / overlap;
289  weight = csi * eta;
290  interpolation *= weight;
291  dval += interpolation;
292  }
293  else if ((Y >= Overlap.S) &&
294  (Y <= Overlap.N)) { /* (2) */
295  weight = (X - General.W) / overlap;
296  interpolation *= weight;
297  dval += interpolation;
298  }
299  }
300  else if ((X >= Overlap.W) && (X <= Overlap.E)) {
301  if ((Y > Overlap.N) && (Y < General.N)) { /* (3) */
302  weight = (General.N - Y) / overlap;
303  interpolation *= weight;
304  dval += interpolation;
305  }
306  else if ((Y < Overlap.S) && (Y > General.S)) { /* (1) */
307  weight = (Y - General.S) / overlap;
308  interpolation *= weight;
309  dval = interpolation;
310  }
311  }
312  }
313  Segment_put(out_seg, &dval, row, col);
314  }
315  } /* END COL */
316  } /* END ROW */
317  return 1;
318 }
double dataInterpolateBicubic(double x, double y, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, double *parVect)
Definition: InterpSpline.c:497
double dataInterpolateBilin(double x, double y, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, double *parVect)
Definition: InterpSpline.c:604
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:159
void P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration, struct bound_box General, struct bound_box Overlap, double **obs, double *param, int *line_num, double pe, double pn, double overlap, int nsplx, int nsply, int num_points, int bilin, struct line_cats *categories, dbDriver *driver, double mean, char *tab_name)
Definition: lidar/raster.c:8
int P_Regular_Points(struct Cell_head *Elaboration, struct Cell_head *Original, struct bound_box General, struct bound_box Overlap, SEGMENT *out_seg, double *param, double passoN, double passoE, double overlap, double mean, int nsplx, int nsply, int nrows, int ncols, int bilin)
Definition: lidar/raster.c:195
int Segment_get(SEGMENT *SEG, void *buf, off_t row, off_t col)
Get value from segment file.
Definition: segment/get.c:37
int Segment_put(SEGMENT *SEG, const void *buf, off_t row, off_t col)
Definition: segment/put.c:43
if(nSHPType==SHPT_ARCM||nSHPType==SHPT_POINTM||nSHPType==SHPT_POLYGONM||nSHPType==SHPT_MULTIPOINTM)
Definition: shpopen.c:1167
Definition: driver.h:21
#define X(j)
#define Y(j)