GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
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/*------------------------------------------------------------------------------------------------*/
8void 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/*------------------------------------------------------------------------------------------------*/
195int 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)
double dataInterpolateBilin(double x, double y, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, double *parVect)
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)
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
#define X(j)
#define Y(j)