GRASS 8 Programmer's Manual 8.6.0dev(2026)-1d1e47ad9d
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#include <grass/vector.h>
7
8/*------------------------------------------------------------------------------------------------*/
9void P_Sparse_Points(struct Map_info *Out, struct Cell_head *Elaboration,
10 struct bound_box General, struct bound_box Overlap,
11 double **obs, double *param, int *line_num, double pe,
12 double pn, double overlap, int nsplx, int nsply,
13 int num_points, int bilin, struct line_cats *categories,
14 dbDriver *driver, double mean, char *tab_name)
15{
16 int i;
17 char buf[1024];
18 dbString sql;
19
20 double interpolation, csi, eta, weight;
21 struct line_pnts *point;
22
23 point = Vect_new_line_struct();
24
25 db_begin_transaction(driver);
26
27 for (i = 0; i < num_points; i++) {
28
29 if (Vect_point_in_box(obs[i][0], obs[i][1], mean,
30 &General)) { /*Here mean is just for asking if obs
31 point is in box */
32
33 if (bilin)
34 interpolation = dataInterpolateBilin(
35 obs[i][0], obs[i][1], pe, pn, nsplx, nsply,
36 Elaboration->west, Elaboration->south, param);
37 else
38 interpolation = dataInterpolateBicubic(
39 obs[i][0], obs[i][1], pe, pn, nsplx, nsply,
40 Elaboration->west, Elaboration->south, param);
41
42 interpolation += mean;
43 Vect_copy_xyz_to_pnts(point, &obs[i][0], &obs[i][1], &interpolation,
44 1);
45
46 if (Vect_point_in_box(obs[i][0], obs[i][1], interpolation,
47 &Overlap)) { /*(5) */
48 Vect_write_line(Out, GV_POINT, point, categories);
49 }
50 else {
51 db_init_string(&sql);
52
53 snprintf(buf, sizeof(buf), "INSERT INTO %s (ID, X, Y, Interp)",
54 tab_name);
55 db_append_string(&sql, buf);
56
57 snprintf(buf, sizeof(buf), " VALUES (");
58 db_append_string(&sql, buf);
59 snprintf(buf, sizeof(buf), "%d, %f, %f, ", line_num[i],
60 obs[i][0], obs[i][1]);
61 db_append_string(&sql, buf);
62
63 if ((*point->x > Overlap.E) && (*point->x < General.E)) {
64 if ((*point->y > Overlap.N) &&
65 (*point->y < General.N)) { /*(3) */
66 csi = (General.E - *point->x) / overlap;
67 eta = (General.N - *point->y) / overlap;
68 weight = csi * eta;
69 *point->z = weight * interpolation;
70
71 snprintf(buf, sizeof(buf), "%lf", *point->z);
72 db_append_string(&sql, buf);
73 snprintf(buf, sizeof(buf), ")");
74 db_append_string(&sql, buf);
75
76 if (db_execute_immediate(driver, &sql) != DB_OK)
77 G_fatal_error(_("Unable to access table <%s>"),
78 tab_name);
79 }
80 else if ((*point->y < Overlap.S) &&
81 (*point->y > General.S)) { /*(1) */
82 csi = (General.E - *point->x) / overlap;
83 eta = (*point->y - General.S) / overlap;
84 weight = csi * eta;
85 *point->z = weight * interpolation;
86
87 snprintf(buf, sizeof(buf), "%lf", *point->z);
88 db_append_string(&sql, buf);
89 snprintf(buf, sizeof(buf), ")");
90 db_append_string(&sql, buf);
91
92 if (db_execute_immediate(driver, &sql) != DB_OK)
93 G_fatal_error(_("Unable to access table <%s>"),
94 tab_name);
95 }
96 else if ((*point->y <= Overlap.N) &&
97 (*point->y >= Overlap.S)) { /*(1) */
98 weight = (General.E - *point->x) / overlap;
99 *point->z = weight * interpolation;
100
101 snprintf(buf, sizeof(buf), "%lf", *point->z);
102 db_append_string(&sql, buf);
103 snprintf(buf, sizeof(buf), ")");
104 db_append_string(&sql, buf);
105
106 if (db_execute_immediate(driver, &sql) != DB_OK)
107 G_fatal_error(_("Unable to access table <%s>"),
108 tab_name);
109 }
110 }
111 else if ((*point->x < Overlap.W) && (*point->x > General.W)) {
112 if ((*point->y > Overlap.N) &&
113 (*point->y < General.N)) { /*(4) */
114 csi = (*point->x - General.W) / overlap;
115 eta = (General.N - *point->y) / overlap;
116 weight = eta * csi;
117 *point->z = weight * interpolation;
118
119 snprintf(buf, sizeof(buf), "%lf", *point->z);
120 db_append_string(&sql, buf);
121 snprintf(buf, sizeof(buf), ")");
122 db_append_string(&sql, buf);
123
124 if (db_execute_immediate(driver, &sql) != DB_OK)
125 G_fatal_error(_("Unable to access table <%s>"),
126 tab_name);
127 }
128 else if ((*point->y < Overlap.S) &&
129 (*point->y > General.S)) { /*(2) */
130 csi = (*point->x - General.W) / overlap;
131 eta = (*point->y - General.S) / overlap;
132 weight = csi * eta;
133 *point->z = weight * interpolation;
134
135 snprintf(buf, sizeof(buf), "%lf", *point->z);
136 db_append_string(&sql, buf);
137 snprintf(buf, sizeof(buf), ")");
138 db_append_string(&sql, buf);
139
140 if (db_execute_immediate(driver, &sql) != DB_OK)
141 G_fatal_error(_("Unable to access table <%s>"),
142 tab_name);
143 }
144 else if ((*point->y >= Overlap.S) &&
145 (*point->y <= Overlap.N)) { /*(2) */
146 weight = (*point->x - General.W) / overlap;
147 *point->z = weight * interpolation;
148
149 snprintf(buf, sizeof(buf), "%lf", *point->z);
150 db_append_string(&sql, buf);
151 snprintf(buf, sizeof(buf), ")");
152 db_append_string(&sql, buf);
153
154 if (db_execute_immediate(driver, &sql) != DB_OK)
155 G_fatal_error(_("Unable to access table <%s>"),
156 tab_name);
157 }
158 }
159 else if ((*point->x >= Overlap.W) && (*point->x <= Overlap.E)) {
160 if ((*point->y > Overlap.N) &&
161 (*point->y < General.N)) { /*(3) */
162 weight = (General.N - *point->y) / overlap;
163 *point->z = weight * interpolation;
164
165 snprintf(buf, sizeof(buf), "%lf", *point->z);
166 db_append_string(&sql, buf);
167 snprintf(buf, sizeof(buf), ")");
168 db_append_string(&sql, buf);
169
170 if (db_execute_immediate(driver, &sql) != DB_OK)
171 G_fatal_error(_("Unable to access table <%s>"),
172 tab_name);
173 }
174 else if ((*point->y < Overlap.S) &&
175 (*point->y > General.S)) { /*(1) */
176 weight = (*point->y - General.S) / overlap;
177 *point->z = (1 - weight) * interpolation;
178
179 snprintf(buf, sizeof(buf), "%lf", *point->z);
180 db_append_string(&sql, buf);
181 snprintf(buf, sizeof(buf), ")");
182 db_append_string(&sql, buf);
183
184 if (db_execute_immediate(driver, &sql) != DB_OK)
185 G_fatal_error(_("Unable to access table <%s>"),
186 tab_name);
187 }
188 }
189 }
190 }
191 /*IF*/ }
192 /*FOR*/ db_commit_transaction(driver);
193 Vect_destroy_line_struct(point);
194
195 return;
196}
197
198/*------------------------------------------------------------------------------------------------*/
199int P_Regular_Points(struct Cell_head *Elaboration, struct Cell_head *Original,
200 struct bound_box General, struct bound_box Overlap,
201 SEGMENT *out_seg, double *param, double passoN,
202 double passoE, double overlap, double mean, int nsplx,
203 int nsply, int nrows, int ncols, int bilin)
204{
205
206 int col, row, startcol, endcol, startrow, endrow;
207 double X, Y, interpolation, weight, csi, eta, dval;
208
209 /* G_get_window(&Original); */
210 if (Original->north > General.N)
211 startrow = (Original->north - General.N) / Original->ns_res - 1;
212 else
213 startrow = 0;
214 if (Original->north > General.S) {
215 endrow = (Original->north - General.S) / Original->ns_res + 1;
216 if (endrow > nrows)
217 endrow = nrows;
218 }
219 else
220 endrow = nrows;
221 if (General.W > Original->west)
222 startcol = (General.W - Original->west) / Original->ew_res - 1;
223 else
224 startcol = 0;
225 if (General.E > Original->west) {
226 endcol = (General.E - Original->west) / Original->ew_res + 1;
227 if (endcol > ncols)
228 endcol = ncols;
229 }
230 else
231 endcol = ncols;
232
233 for (row = startrow; row < endrow; row++) {
234 for (col = startcol; col < endcol; col++) {
235
236 X = Rast_col_to_easting((double)(col) + 0.5, Original);
237 Y = Rast_row_to_northing((double)(row) + 0.5, Original);
238
239 if (Vect_point_in_box(X, Y, mean,
240 &General)) { /* Here, mean is just for asking
241 if obs point is in box */
242
243 if (bilin)
244 interpolation = dataInterpolateBilin(
245 X, Y, passoE, passoN, nsplx, nsply, Elaboration->west,
246 Elaboration->south, param);
247 else
248 interpolation = dataInterpolateBicubic(
249 X, Y, passoE, passoN, nsplx, nsply, Elaboration->west,
250 Elaboration->south, param);
251
252 interpolation += mean;
253
254 if (Vect_point_in_box(X, Y, interpolation,
255 &Overlap)) { /* (5) */
256 dval = interpolation;
257 }
258 else {
259 Segment_get(out_seg, &dval, row, col);
260 if ((X > Overlap.E) && (X < General.E)) {
261 if ((Y > Overlap.N) && (Y < General.N)) { /* (3) */
262 csi = (General.E - X) / overlap;
263 eta = (General.N - Y) / overlap;
264 weight = csi * eta;
265 interpolation *= weight;
266 dval += interpolation;
267 }
268 else if ((Y < Overlap.S) && (Y > General.S)) { /* (1) */
269 csi = (General.E - X) / overlap;
270 eta = (Y - General.S) / overlap;
271 weight = csi * eta;
272 interpolation *= weight;
273 dval = interpolation;
274 }
275 else if ((Y >= Overlap.S) &&
276 (Y <= Overlap.N)) { /* (1) */
277 weight = (General.E - X) / overlap;
278 interpolation *= weight;
279 dval = interpolation;
280 }
281 }
282 else if ((X < Overlap.W) && (X > General.W)) {
283 if ((Y > Overlap.N) && (Y < General.N)) { /* (4) */
284 csi = (X - General.W) / overlap;
285 eta = (General.N - Y) / overlap;
286 weight = eta * csi;
287 interpolation *= weight;
288 dval += interpolation;
289 }
290 else if ((Y < Overlap.S) && (Y > General.S)) { /* (2) */
291 csi = (X - General.W) / overlap;
292 eta = (Y - General.S) / overlap;
293 weight = csi * eta;
294 interpolation *= weight;
295 dval += interpolation;
296 }
297 else if ((Y >= Overlap.S) &&
298 (Y <= Overlap.N)) { /* (2) */
299 weight = (X - General.W) / overlap;
300 interpolation *= weight;
301 dval += interpolation;
302 }
303 }
304 else if ((X >= Overlap.W) && (X <= Overlap.E)) {
305 if ((Y > Overlap.N) && (Y < General.N)) { /* (3) */
306 weight = (General.N - Y) / overlap;
307 interpolation *= weight;
308 dval += interpolation;
309 }
310 else if ((Y < Overlap.S) && (Y > General.S)) { /* (1) */
311 weight = (Y - General.S) / overlap;
312 interpolation *= weight;
313 dval = interpolation;
314 }
315 }
316 }
317 Segment_put(out_seg, &dval, row, col);
318 }
319 } /* END COL */
320 } /* END ROW */
321 return 1;
322}
dataInterpolateBicubic
double dataInterpolateBicubic(double x, double y, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, double *parVect)
Definition InterpSpline.c:497
dataInterpolateBilin
double dataInterpolateBilin(double x, double y, double deltaX, double deltaY, int xNum, int yNum, double xMin, double yMin, double *parVect)
Definition InterpSpline.c:604
AMI_STREAM
Definition ami_stream.h:153
DB_OK
#define DB_OK
Definition dbmi.h:71
db_commit_transaction
int db_commit_transaction(dbDriver *)
Commit transaction.
Definition c_execute.c:82
db_begin_transaction
int db_begin_transaction(dbDriver *)
Begin transaction.
Definition c_execute.c:56
db_execute_immediate
int db_execute_immediate(dbDriver *, dbString *)
Execute SQL statements.
Definition c_execute.c:27
db_init_string
void db_init_string(dbString *)
Initialize dbString.
Definition string.c:25
db_append_string
int db_append_string(dbString *, const char *)
Append string to dbString.
Definition string.c:205
G_fatal_error
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
Rast_row_to_northing
double Rast_row_to_northing(double, const struct Cell_head *)
Row to northing.
Definition raster/window.c:244
Rast_col_to_easting
double Rast_col_to_easting(double, const struct Cell_head *)
Column to easting.
Definition raster/window.c:264
Segment_get
int Segment_get(SEGMENT *, void *, off_t, off_t)
Get value from segment file.
Definition segment/get.c:36
Segment_put
int Segment_put(SEGMENT *, const void *, off_t, off_t)
Definition segment/put.c:42
Vect_destroy_line_struct
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
Definition line.c:77
Vect_copy_xyz_to_pnts
int Vect_copy_xyz_to_pnts(struct line_pnts *, const double *, const double *, const double *, int)
Copy points from array to line_pnts structure.
Definition line.c:99
Vect_write_line
off_t Vect_write_line(struct Map_info *, int, const struct line_pnts *, const struct line_cats *)
Writes a new feature.
Definition vector/Vlib/write.c:168
Vect_point_in_box
int Vect_point_in_box(double, double, double, const struct bound_box *)
Tests if point is in 3D box.
Definition vector/Vlib/box.c:48
Vect_new_line_struct
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition line.c:45
GV_POINT
#define GV_POINT
Feature types used in memory on run time (may change)
Definition dig_defines.h:183
_
#define _(str)
Definition glocale.h:10
mean
float mean(IClass_statistics *statistics, int band)
Helper function for computing mean.
Definition iclass_statistics.c:369
P_Sparse_Points
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:9
P_Regular_Points
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:199
X
#define X
Definition ogsf.h:140
Y
#define Y
Definition ogsf.h:141
Cell_head
2D/3D raster map header (used also for region)
Definition gis.h:446
Map_info
Vector map info.
Definition dig_structs.h:1231
_db_driver
Definition dbmi.h:165
_db_string
Definition dbmi.h:147
bound_box
Bounding box.
Definition dig_structs.h:62
line_cats
Feature category info.
Definition dig_structs.h:1665
line_pnts
Feature geometry info - coordinates.
Definition dig_structs.h:1639
line_pnts::y
double * y
Array of Y coordinates.
Definition dig_structs.h:1647
line_pnts::x
double * x
Array of X coordinates.
Definition dig_structs.h:1643
line_pnts::z
double * z
Array of Z coordinates.
Definition dig_structs.h:1651