1 #include <grass/config.h>
6 #include <grass/lidar.h>
54 P_set_regions(
struct Cell_head *Elaboration,
struct bound_box * General,
55 struct bound_box * Overlap,
struct Reg_dimens dim,
int type)
59 struct Cell_head orig;
67 Elaboration->south = Elaboration->north - dim.
sn_size;
68 General->N = Elaboration->north - dim.
edge_h;
69 General->S = Elaboration->south + dim.
edge_h;
70 Overlap->N = General->N - dim.
overlap;
71 Overlap->S = General->S + dim.
overlap;
77 Elaboration->east = Elaboration->west + dim.
ew_size;
78 General->W = Elaboration->west + dim.
edge_v;
79 General->E = Elaboration->east - dim.
edge_v;
80 Overlap->W = General->W + dim.
overlap;
81 Overlap->E = General->E - dim.
overlap;
85 Elaboration->north = orig.north + 2 * dim.
edge_h;
86 Elaboration->south = Elaboration->north - dim.
sn_size;
87 General->N = orig.north;
88 General->S = Elaboration->south + dim.
edge_h;
89 Overlap->N = General->N;
90 Overlap->S = General->S + dim.
overlap;
94 Elaboration->south = orig.south - 2 * dim.
edge_h;
95 General->S = orig.south;
96 Overlap->S = General->S;
100 Elaboration->west = orig.west - 2 * dim.
edge_v;
101 Elaboration->east = Elaboration->west + dim.
ew_size;
102 General->W = orig.west;
103 General->E = Elaboration->east - dim.
edge_v;
104 Overlap->W = General->W;
105 Overlap->E = General->E - dim.
overlap;
109 Elaboration->east = orig.east + 2 * dim.
edge_v;
110 General->E = orig.east;
111 Overlap->E = General->E;
121 int total_splines, edge_splines, n_windows;
122 int lastsplines, lastsplines_min, lastsplines_max;
123 double E_extension, N_extension, edgeE, edgeN;
124 struct Cell_head orig;
129 E_extension = orig.east - orig.west;
130 N_extension = orig.north - orig.south;
139 total_splines = ceil(E_extension / pe);
140 edge_splines = edgeE / pe;
141 n_windows = E_extension / edgeE;
147 lastsplines = total_splines - edge_splines * n_windows;
148 while (lastsplines > lastsplines_max || lastsplines < lastsplines_min) {
153 edge_splines = edgeE / pe;
154 n_windows = E_extension / edgeE;
155 lastsplines = total_splines - edge_splines * n_windows;
161 total_splines = ceil(N_extension / pn);
162 edge_splines = edgeN / pn;
163 n_windows = N_extension / edgeN;
169 lastsplines = total_splines - edge_splines * n_windows;
170 while (lastsplines > lastsplines_max || lastsplines < lastsplines_min) {
175 edge_splines = edgeN / pn;
176 n_windows = N_extension / edgeN;
177 lastsplines = total_splines - edge_splines * n_windows;
213 return (2 * nsplines + 1);
216 return (4 * nsplines + 3);
224 int i, mean_count = 0;
226 struct bound_box mean_box;
228 Vect_region_box(Elaboration, &mean_box);
234 for (i = 0; i < npoints; i++) {
235 if (Vect_point_in_box
236 (obs[i].coordX, obs[i].coordY, obs[i].coordZ, &mean_box)) {
244 mean /= (double)mean_count;
251 int type, npoints = 0;
252 double xmin = 0, xmax = 0, ymin = 0, ymax = 0;
254 struct line_pnts *points;
255 struct line_cats *categories;
256 struct bound_box region_box;
257 struct Cell_head orig;
260 Vect_region_box(&orig, ®ion_box);
262 points = Vect_new_line_struct();
263 categories = Vect_new_cats_struct();
266 while ((type = Vect_read_next_line(Map, points, categories)) > 0) {
267 if (!(type & GV_POINT))
272 if (points->z !=
NULL)
278 if (Vect_point_in_box(x, y, z, ®ion_box)) {
299 *dist = sqrt(((xmax - xmin) * (ymax - ymin)) / npoints);
301 *dens = npoints / ((xmax - xmin) * (ymax - ymin));
310 struct Cell_head *Elaboration,
311 int *num_points,
int dim_vect,
314 int line_num, pippo, npoints,
cat, type;
317 struct line_pnts *points;
318 struct line_cats *categories;
319 struct bound_box elaboration_box;
322 obs = (
struct Point *)G_calloc(pippo,
sizeof(
struct Point));
324 points = Vect_new_line_struct();
325 categories = Vect_new_cats_struct();
328 Vect_region_box(Elaboration, &elaboration_box);
334 while ((type = Vect_read_next_line(Map, points, categories)) > 0) {
336 if (!(type & GV_POINT))
343 if (points->z !=
NULL)
349 if (Vect_point_in_box(x, y, z, &elaboration_box)) {
351 if (npoints >= pippo) {
354 (
struct Point *)G_realloc((
void *)obs,
356 sizeof(
struct Point));
360 obs[npoints - 1].
coordX = x;
361 obs[npoints - 1].
coordY = y;
362 obs[npoints - 1].
coordZ = z;
363 obs[npoints - 1].
lineID = line_num;
365 Vect_cat_get(categories, layer, &cat);
366 obs[npoints - 1].
cat =
cat;
369 Vect_destroy_line_struct(points);
370 Vect_destroy_cats_struct(categories);
372 *num_points = npoints;
377 struct Cell_head *Elaboration,
378 struct Cell_head *Original,
379 int *num_points,
int dim_vect)
381 int col, row, startcol, endcol, startrow, endrow, nrows, ncols;
385 struct bound_box elaboration_box;
388 obs = (
struct Point *)G_calloc(pippo,
sizeof(
struct Point));
391 Vect_region_box(Elaboration, &elaboration_box);
394 nrows = Original->rows;
395 ncols = Original->cols;
397 if (Original->north > Elaboration->north)
398 startrow = (Original->north - Elaboration->north) / Original->ns_res - 1;
401 if (Original->north > Elaboration->south) {
402 endrow = (Original->north - Elaboration->south) / Original->ns_res + 1;
408 if (Elaboration->west > Original->west)
409 startcol = (Elaboration->west - Original->west) / Original->ew_res - 1;
412 if (Elaboration->east > Original->west) {
413 endcol = (Elaboration->east - Original->west) / Original->ew_res + 1;
420 for (row = startrow; row < endrow; row++) {
421 for (col = startcol; col < endcol; col++) {
425 if (!Rast_is_d_null_value(&z)) {
427 x = Rast_col_to_easting((
double)(col) + 0.5, Original);
428 y = Rast_row_to_northing((
double)(row) + 0.5, Original);
430 if (Vect_point_in_box(x, y, 0, &elaboration_box)) {
432 if (npoints >= pippo) {
435 (
struct Point *)G_realloc((
void *)obs,
437 sizeof(
struct Point));
441 obs[npoints - 1].
coordX = x;
442 obs[npoints - 1].
coordY = y;
443 obs[npoints - 1].
coordZ = z;
449 *num_points = npoints;
456 dbTable *auxiliar_tab;
459 auxiliar_tab = db_alloc_table(2);
460 db_set_table_name(auxiliar_tab, tab_name);
461 db_set_table_description(auxiliar_tab,
462 "Intermediate interpolated values");
464 column = db_get_table_column(auxiliar_tab, 0);
465 db_set_column_name(column,
"ID");
466 db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);
468 column = db_get_table_column(auxiliar_tab, 1);
469 db_set_column_name(column,
"Interp");
470 db_set_column_sqltype(column, DB_SQL_TYPE_REAL);
472 if (db_create_table(driver, auxiliar_tab) == DB_OK) {
473 G_debug(1, _(
"<%s> created in database."), tab_name);
477 G_warning(_(
"<%s> has not been created in database."), tab_name);
485 dbTable *auxiliar_tab;
488 auxiliar_tab = db_alloc_table(4);
489 db_set_table_name(auxiliar_tab, tab_name);
490 db_set_table_description(auxiliar_tab,
491 "Intermediate interpolated values");
493 column = db_get_table_column(auxiliar_tab, 0);
494 db_set_column_name(column,
"ID");
495 db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);
497 column = db_get_table_column(auxiliar_tab, 1);
498 db_set_column_name(column,
"Interp");
499 db_set_column_sqltype(column, DB_SQL_TYPE_REAL);
501 column = db_get_table_column(auxiliar_tab, 2);
502 db_set_column_name(column,
"X");
503 db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
505 column = db_get_table_column(auxiliar_tab, 3);
506 db_set_column_name(column,
"Y");
507 db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
509 if (db_create_table(driver, auxiliar_tab) == DB_OK) {
510 G_debug(1, _(
"<%s> created in database."), tab_name);
514 G_warning(_(
"<%s> has not been created in database."), tab_name);
524 db_init_string(&drop);
525 db_append_string(&drop,
"drop table ");
526 db_append_string(&drop, tab_name);
527 return db_execute_immediate(driver, &drop);
533 int ncols, col, nrows, row;
536 nrows = Rast_window_rows();
537 ncols = Rast_window_cols();
539 raster = Rast_allocate_buf(DCELL_TYPE);
541 for (row = 0; row < nrows; row++) {
544 Rast_set_d_null_value(raster, ncols);
546 for (col = 0, ptr = raster; col < ncols;
548 Rast_set_d_value(ptr, (DCELL) (matrix[row][col]), DCELL_TYPE);
550 Rast_put_d_row(fd, raster);
561 int more, line_num, type,
count = 0;
564 struct line_pnts *point;
565 struct line_cats *cat;
574 point = Vect_new_line_struct();
575 cat = Vect_new_cats_struct();
577 db_init_string(&sql);
578 db_zero_string(&sql);
580 sprintf(buf,
"select ID, X, Y, sum(Interp) from %s group by ID, X, Y",
583 db_append_string(&sql, buf);
584 db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL);
586 while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
588 table = db_get_cursor_table(&cursor);
590 column = db_get_table_column(table, 0);
591 type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
592 if (type == DB_C_TYPE_INT)
593 value = db_get_column_value(column);
596 line_num = db_get_value_int(value);
598 column = db_get_table_column(table, 1);
599 type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
600 if (type == DB_C_TYPE_DOUBLE)
601 value = db_get_column_value(column);
604 coordZ = db_get_value_double(value);
606 column = db_get_table_column(table, 2);
607 type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
608 if (type == DB_C_TYPE_DOUBLE)
609 value = db_get_column_value(column);
612 coordX = db_get_value_double(value);
614 column = db_get_table_column(table, 3);
615 type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
616 if (type == DB_C_TYPE_DOUBLE)
617 value = db_get_column_value(column);
620 coordY = db_get_value_double(value);
622 Vect_copy_xyz_to_pnts(point, &coordX, &coordY, &coordZ, 1);
623 Vect_reset_cats(cat);
624 Vect_cat_set(cat, 1, 1);
625 Vect_write_line(Out, GV_POINT, point, cat);
void G_get_window(struct Cell_head *window)
Get the current region.
int Segment_get(SEGMENT *SEG, void *buf, off_t row, off_t col)
Get value from segment file.
void P_zero_dim(struct Reg_dimens *dim)
struct Point * P_Read_Vector_Region_Map(struct Map_info *Map, struct Cell_head *Elaboration, int *num_points, int dim_vect, int layer)
void P_Aux_to_Vector(struct Map_info *Map, struct Map_info *Out, dbDriver *driver, char *tab_name)
int P_get_edge(int interpolator, struct Reg_dimens *dim, double pe, double pn)
int P_Create_Aux4_Table(dbDriver *driver, char *tab_name)
int G_debug(int level, const char *msg,...)
Print debugging message.
void P_Aux_to_Raster(double **matrix, int fd)
double P_estimate_splinestep(struct Map_info *Map, double *dens, double *dist)
int P_Drop_Aux_Table(dbDriver *driver, char *tab_name)
if(!DBFLoadRecord(psDBF, hEntity)) return NULL
void G_percent(long n, long d, int s)
Print percent complete messages.
int P_Create_Aux2_Table(dbDriver *driver, char *tab_name)
int P_set_dim(struct Reg_dimens *dim, double pe, double pn, int *nsplx, int *nsply)
int P_get_BandWidth(int interpolator, int nsplines)
void G_get_set_window(struct Cell_head *window)
Get the current working window (region)
double P_Mean_Calc(struct Cell_head *Elaboration, struct Point *obs, int npoints)
int P_set_regions(struct Cell_head *Elaboration, struct bound_box *General, struct bound_box *Overlap, struct Reg_dimens dim, int type)
void * G_incr_void_ptr(const void *ptr, size_t size)
Advance void pointer.
void G_warning(const char *msg,...)
Print a warning message to stderr.
struct Point * P_Read_Raster_Region_Map(SEGMENT *in_seg, struct Cell_head *Elaboration, struct Cell_head *Original, int *num_points, int dim_vect)