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)) {
297 Vect_destroy_cats_struct(categories);
298 Vect_destroy_line_struct(points);
302 *dist = sqrt(((xmax - xmin) * (ymax - ymin)) / npoints);
304 *dens = npoints / ((xmax - xmin) * (ymax - ymin));
313 struct Cell_head *Elaboration,
314 int *num_points,
int dim_vect,
317 int line_num, pippo, npoints,
cat, type;
320 struct line_pnts *points;
321 struct line_cats *categories;
322 struct bound_box elaboration_box;
325 obs = (
struct Point *)G_calloc(pippo,
sizeof(
struct Point));
327 points = Vect_new_line_struct();
328 categories = Vect_new_cats_struct();
331 Vect_region_box(Elaboration, &elaboration_box);
337 while ((type = Vect_read_next_line(Map, points, categories)) > 0) {
339 if (!(type & GV_POINT))
346 if (points->z !=
NULL)
352 if (Vect_point_in_box(x, y, z, &elaboration_box)) {
354 if (npoints >= pippo) {
357 (
struct Point *)G_realloc((
void *)obs,
359 sizeof(
struct Point));
364 obs[npoints - 1].
coordY = y;
365 obs[npoints - 1].
coordZ = z;
366 obs[npoints - 1].
lineID = line_num;
368 Vect_cat_get(categories, layer, &cat);
369 obs[npoints - 1].
cat =
cat;
372 Vect_destroy_line_struct(points);
373 Vect_destroy_cats_struct(categories);
375 *num_points = npoints;
380 struct Cell_head *Elaboration,
381 struct Cell_head *Original,
382 int *num_points,
int dim_vect)
384 int col, row, startcol, endcol, startrow, endrow, nrows, ncols;
388 struct bound_box elaboration_box;
391 obs = (
struct Point *)G_calloc(pippo,
sizeof(
struct Point));
394 Vect_region_box(Elaboration, &elaboration_box);
397 nrows = Original->rows;
398 ncols = Original->cols;
400 if (Original->north > Elaboration->north)
401 startrow = (Original->north - Elaboration->north) / Original->ns_res - 1;
404 if (Original->north > Elaboration->south) {
405 endrow = (Original->north - Elaboration->south) / Original->ns_res + 1;
411 if (Elaboration->west > Original->west)
412 startcol = (Elaboration->west - Original->west) / Original->ew_res - 1;
415 if (Elaboration->east > Original->west) {
416 endcol = (Elaboration->east - Original->west) / Original->ew_res + 1;
423 for (row = startrow; row < endrow; row++) {
424 for (col = startcol; col < endcol; col++) {
428 if (!Rast_is_d_null_value(&z)) {
430 x = Rast_col_to_easting((
double)(col) + 0.5, Original);
431 y = Rast_row_to_northing((
double)(row) + 0.5, Original);
433 if (Vect_point_in_box(x, y, 0, &elaboration_box)) {
435 if (npoints >= pippo) {
438 (
struct Point *)G_realloc((
void *)obs,
440 sizeof(
struct Point));
445 obs[npoints - 1].
coordY = y;
446 obs[npoints - 1].
coordZ = z;
452 *num_points = npoints;
459 dbTable *auxiliar_tab;
462 auxiliar_tab = db_alloc_table(2);
463 db_set_table_name(auxiliar_tab, tab_name);
464 db_set_table_description(auxiliar_tab,
465 "Intermediate interpolated values");
467 column = db_get_table_column(auxiliar_tab, 0);
468 db_set_column_name(column,
"ID");
469 db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);
471 column = db_get_table_column(auxiliar_tab, 1);
472 db_set_column_name(column,
"Interp");
473 db_set_column_sqltype(column, DB_SQL_TYPE_REAL);
475 if (db_create_table(driver, auxiliar_tab) == DB_OK) {
476 G_debug(1, _(
"<%s> created in database."), tab_name);
480 G_warning(_(
"<%s> has not been created in database."), tab_name);
488 dbTable *auxiliar_tab;
491 auxiliar_tab = db_alloc_table(4);
492 db_set_table_name(auxiliar_tab, tab_name);
493 db_set_table_description(auxiliar_tab,
494 "Intermediate interpolated values");
496 column = db_get_table_column(auxiliar_tab, 0);
497 db_set_column_name(column,
"ID");
498 db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);
500 column = db_get_table_column(auxiliar_tab, 1);
501 db_set_column_name(column,
"Interp");
502 db_set_column_sqltype(column, DB_SQL_TYPE_REAL);
504 column = db_get_table_column(auxiliar_tab, 2);
505 db_set_column_name(column,
"X");
506 db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
508 column = db_get_table_column(auxiliar_tab, 3);
509 db_set_column_name(column,
"Y");
510 db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
512 if (db_create_table(driver, auxiliar_tab) == DB_OK) {
513 G_debug(1, _(
"<%s> created in database."), tab_name);
517 G_warning(_(
"<%s> has not been created in database."), tab_name);
527 db_init_string(&drop);
528 db_append_string(&drop,
"drop table ");
529 db_append_string(&drop, tab_name);
530 return db_execute_immediate(driver, &drop);
536 int ncols, col, nrows, row;
539 nrows = Rast_window_rows();
540 ncols = Rast_window_cols();
542 raster = Rast_allocate_buf(DCELL_TYPE);
544 for (row = 0; row < nrows; row++) {
547 Rast_set_d_null_value(raster, ncols);
549 for (col = 0, ptr = raster; col < ncols;
551 Rast_set_d_value(ptr, (DCELL) (matrix[row][col]), DCELL_TYPE);
553 Rast_put_d_row(fd, raster);
564 int more, line_num, type,
count = 0;
567 struct line_pnts *point;
568 struct line_cats *cat;
577 point = Vect_new_line_struct();
578 cat = Vect_new_cats_struct();
580 db_init_string(&sql);
581 db_zero_string(&sql);
583 sprintf(buf,
"select ID, X, Y, sum(Interp) from %s group by ID, X, Y",
586 db_append_string(&sql, buf);
587 db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL);
589 while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
591 table = db_get_cursor_table(&cursor);
593 column = db_get_table_column(table, 0);
594 type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
595 if (type == DB_C_TYPE_INT)
596 value = db_get_column_value(column);
599 line_num = db_get_value_int(value);
601 column = db_get_table_column(table, 1);
602 type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
603 if (type == DB_C_TYPE_DOUBLE)
604 value = db_get_column_value(column);
607 coordZ = db_get_value_double(value);
609 column = db_get_table_column(table, 2);
610 type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
611 if (type == DB_C_TYPE_DOUBLE)
612 value = db_get_column_value(column);
615 coordX = db_get_value_double(value);
617 column = db_get_table_column(table, 3);
618 type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
619 if (type == DB_C_TYPE_DOUBLE)
620 value = db_get_column_value(column);
623 coordY = db_get_value_double(value);
625 Vect_copy_xyz_to_pnts(point, &coordX, &coordY, &coordZ, 1);
626 Vect_reset_cats(cat);
627 Vect_cat_set(cat, 1, 1);
628 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)