Loading...
Searching...
No Matches
zones.c
Go to the documentation of this file.
#include <grass/config.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <grass/lidar.h>
/*----------------------------------------------------------------------------------------*/
void P_zero_dim(struct Reg_dimens *dim)
{
dim->edge_h = 0.0;
dim->edge_v = 0.0;
dim->overlap = 0.0;
dim->sn_size = 0.0;
dim->ew_size = 0.0;
return;
}
/*----------------------------------------------------------------------------------------*/
/*
--------------------------------------------
| Elaboration region |
| ------------------------------------ |
| | General region | |
| | ---------------------------- | |
| | | | | |
| | | | | |
| | | | | |
| | | Overlap region | | |
| | | | | |
| | | | | |
| | | | | |
| | ---------------------------- | |
| | | |
| ------------------------------------ |
| |
--------------------------------------------
The terminology is misleading:
The Overlap region does NOT overlap with neighbouring segments,
but the Elaboration and General region do overlap
Elaboration is used for interpolation
Interpolated points in Elaboration but outside General are discarded
Interpolated points in General but outside Overlap are weighed by
their distance to Overlap and summed up
Interpolated points in Overlap are taken as they are
The buffer zones Elaboration - General and General - Overlap must be
large enough to avoid artifacts
*/
int P_set_regions(struct Cell_head *Elaboration, struct bound_box *General,
struct bound_box *Overlap, struct Reg_dimens dim, int type)
{
/* Set the Elaboration, General, and Overlap region limits
* Returns 0 on success; -1 on failure*/
struct Cell_head orig;
G_get_window(&orig);
switch (type) {
case GENERAL_ROW: /* General case N-S direction */
Elaboration->north =
Elaboration->south + dim.overlap + (2 * dim.edge_h);
Elaboration->south = Elaboration->north - dim.sn_size;
General->N = Elaboration->north - dim.edge_h;
General->S = Elaboration->south + dim.edge_h;
Overlap->N = General->N - dim.overlap;
Overlap->S = General->S + dim.overlap;
return 0;
case GENERAL_COLUMN: /* General case E-W direction */
Elaboration->west = Elaboration->east - dim.overlap - (2 * dim.edge_v);
Elaboration->east = Elaboration->west + dim.ew_size;
General->W = Elaboration->west + dim.edge_v;
General->E = Elaboration->east - dim.edge_v;
Overlap->W = General->W + dim.overlap;
Overlap->E = General->E - dim.overlap;
return 0;
case FIRST_ROW: /* Just started with first row */
Elaboration->north = orig.north + 2 * dim.edge_h;
Elaboration->south = Elaboration->north - dim.sn_size;
General->N = orig.north;
General->S = Elaboration->south + dim.edge_h;
Overlap->N = General->N;
Overlap->S = General->S + dim.overlap;
return 0;
case LAST_ROW: /* Reached last row */
Elaboration->south = orig.south - 2 * dim.edge_h;
General->S = orig.south;
Overlap->S = General->S;
return 0;
case FIRST_COLUMN: /* Just started with first column */
Elaboration->west = orig.west - 2 * dim.edge_v;
Elaboration->east = Elaboration->west + dim.ew_size;
General->W = orig.west;
General->E = Elaboration->east - dim.edge_v;
Overlap->W = General->W;
Overlap->E = General->E - dim.overlap;
return 0;
case LAST_COLUMN: /* Reached last column */
Elaboration->east = orig.east + 2 * dim.edge_v;
General->E = orig.east;
Overlap->E = General->E;
return 0;
}
return -1;
}
/*----------------------------------------------------------------------------------------*/
int P_set_dim(struct Reg_dimens *dim, double pe, double pn, int *nsplx,
int *nsply)
{
int total_splines, edge_splines, n_windows;
int lastsplines, lastsplines_min, lastsplines_max;
double E_extension, N_extension, edgeE, edgeN;
struct Cell_head orig;
int ret = 0;
G_get_window(&orig);
E_extension = orig.east - orig.west;
N_extension = orig.north - orig.south;
dim->ew_size = *nsplx * pe;
dim->sn_size = *nsply * pn;
edgeE = dim->ew_size - dim->overlap - 2 * dim->edge_v;
edgeN = dim->sn_size - dim->overlap - 2 * dim->edge_h;
/* number of moving windows: E_extension / edgeE */
/* remaining steps: total steps - (floor(E_extension / edgeE) * E_extension)
* / passoE */
/* remaining steps must be larger than edge_v + overlap + half of overlap
* window */
total_splines = ceil(E_extension / pe);
edge_splines = edgeE / pe;
n_windows = E_extension / edgeE; /* without last one */
if (n_windows > 0) {
/* min size of the last overlap window = half of current overlap window
*/
/* max size of the last overlap window = elaboration - 3 * edge -
* overlap */
lastsplines_min =
ceil((dim->ew_size / 2.0 - (dim->edge_v + dim->overlap)) / pe);
lastsplines_max =
ceil((dim->ew_size - 3 * dim->edge_v - dim->overlap) / pe);
lastsplines = total_splines - edge_splines * n_windows;
while (lastsplines > lastsplines_max || lastsplines < lastsplines_min) {
*nsplx -= 1;
dim->ew_size = *nsplx * pe;
edgeE = dim->ew_size - dim->overlap - 2 * dim->edge_v;
edge_splines = edgeE / pe;
n_windows = E_extension / edgeE; /* without last one */
lastsplines = total_splines - edge_splines * n_windows;
if (ret == 0)
ret = 1;
}
}
total_splines = ceil(N_extension / pn);
edge_splines = edgeN / pn;
n_windows = N_extension / edgeN; /* without last one */
if (n_windows > 0) {
/* min size of the last overlap window = half of current overlap window
*/
/* max size of the last overlap window = elaboration - 3 * edge -
* overlap */
lastsplines_min =
ceil((dim->sn_size / 2.0 - (dim->edge_h + dim->overlap)) / pn);
lastsplines_max =
ceil((dim->sn_size - 3 * dim->edge_h - dim->overlap) / pn);
lastsplines = total_splines - edge_splines * n_windows;
while (lastsplines > lastsplines_max || lastsplines < lastsplines_min) {
*nsply -= 1;
dim->sn_size = *nsply * pn;
edgeN = dim->sn_size - dim->overlap - 2 * dim->edge_h;
edge_splines = edgeN / pn;
n_windows = N_extension / edgeN; /* without last one */
lastsplines = total_splines - edge_splines * n_windows;
if (ret < 2)
ret += 2;
}
}
return ret;
}
/*----------------------------------------------------------------------------------------*/
int P_get_edge(int interpolator, struct Reg_dimens *dim, double pe, double pn)
{
/* Set the edge regions dimension
* Returns 1 on success of bilinear; 2 on success of bicubic, 0 on failure
*/
if (interpolator == P_BILINEAR) {
/* in case of edge artifacts, increase as multiples of 3 */
dim->edge_v = 9 * pe;
dim->edge_h = 9 * pn;
return 1;
}
else if (interpolator == P_BICUBIC) {
/* in case of edge artifacts, increase as multiples of 4 */
dim->edge_v = 12 * pe; /*3 */
dim->edge_h = 12 * pn;
return 2;
}
else
return 0; /* The interpolator is neither bilinear nor bicubic!! */
}
/*----------------------------------------------------------------------------------------*/
int P_get_BandWidth(int interpolator, int nsplines)
{
/* Returns the interpolation matrixes BandWidth dimension */
if (interpolator == P_BILINEAR) {
return (2 * nsplines + 1);
}
else {
return (4 * nsplines + 3);
}
}
/*----------------------------------------------------------------------------------------*/
double P_Mean_Calc(struct Cell_head *Elaboration, struct Point *obs,
int npoints)
{
int i, mean_count = 0;
double mean = 0.0;
struct bound_box mean_box;
Vect_region_box(Elaboration, &mean_box);
mean_box.W -= CONTOUR;
mean_box.E += CONTOUR;
mean_box.N += CONTOUR;
mean_box.S -= CONTOUR;
for (i = 0; i < npoints; i++) { /* */
if (Vect_point_in_box(obs[i].coordX, obs[i].coordY, obs[i].coordZ,
&mean_box)) {
mean_count++;
mean += obs[i].coordZ;
}
}
if (mean_count == 0)
mean = .0;
else
mean /= (double)mean_count;
return mean;
}
double P_estimate_splinestep(struct Map_info *Map, double *dens, double *dist)
{
int type, npoints = 0;
double xmin = 0, xmax = 0, ymin = 0, ymax = 0;
double x, y, z;
struct line_pnts *points;
struct line_cats *categories;
struct bound_box region_box;
struct Cell_head orig;
G_get_set_window(&orig);
Vect_region_box(&orig, &region_box);
points = Vect_new_line_struct();
categories = Vect_new_cats_struct();
Vect_rewind(Map);
while ((type = Vect_read_next_line(Map, points, categories)) > 0) {
if (!(type & GV_POINT))
continue;
x = points->x[0];
y = points->y[0];
if (points->z != NULL)
z = points->z[0];
else
z = 0.0;
/* only use points in current region */
if (Vect_point_in_box(x, y, z, &region_box)) {
npoints++;
if (npoints > 1) {
if (xmin > x)
xmin = x;
else if (xmax < x)
xmax = x;
if (ymin > y)
ymin = y;
else if (ymax < y)
ymax = y;
}
else {
xmin = xmax = x;
ymin = ymax = y;
}
}
}
Vect_destroy_cats_struct(categories);
Vect_destroy_line_struct(points);
if (npoints > 0) {
/* estimated average distance between points in map units */
*dist = sqrt(((xmax - xmin) * (ymax - ymin)) / npoints);
/* estimated point density as number of points per square map unit */
*dens = npoints / ((xmax - xmin) * (ymax - ymin));
return 0;
}
else {
return -1;
}
}
struct Point *P_Read_Vector_Region_Map(struct Map_info *Map,
struct Cell_head *Elaboration,
int *num_points, int dim_vect, int layer)
{
int line_num, pippo, npoints, cat, type;
double x, y, z;
struct Point *obs;
struct line_pnts *points;
struct line_cats *categories;
struct bound_box elaboration_box;
pippo = dim_vect;
obs = (struct Point *)G_calloc(pippo, sizeof(struct Point));
points = Vect_new_line_struct();
categories = Vect_new_cats_struct();
/* Reading points inside elaboration zone */
Vect_region_box(Elaboration, &elaboration_box);
npoints = 0;
line_num = 0;
Vect_rewind(Map);
while ((type = Vect_read_next_line(Map, points, categories)) > 0) {
if (!(type & GV_POINT))
continue;
line_num++;
x = points->x[0];
y = points->y[0];
if (points->z != NULL)
z = points->z[0];
else
z = 0.0;
/* Reading and storing points only if in elaboration_reg */
if (Vect_point_in_box(x, y, z, &elaboration_box)) {
npoints++;
if (npoints >= pippo) {
pippo += dim_vect;
obs = (struct Point *)G_realloc(
(void *)obs, (signed int)pippo * sizeof(struct Point));
}
/* Storing observation vector */
obs[npoints - 1].coordX = x;
obs[npoints - 1].coordY = y;
obs[npoints - 1].coordZ = z;
obs[npoints - 1].lineID =
line_num; /* Storing also the line's number */
Vect_cat_get(categories, layer, &cat);
obs[npoints - 1].cat = cat;
}
}
Vect_destroy_line_struct(points);
Vect_destroy_cats_struct(categories);
*num_points = npoints;
return obs;
}
struct Point *P_Read_Raster_Region_Map(SEGMENT *in_seg,
struct Cell_head *Elaboration,
struct Cell_head *Original,
int *num_points, int dim_vect)
{
int col, row, startcol, endcol, startrow, endrow, nrows, ncols;
int pippo, npoints;
double x, y, z;
struct Point *obs;
struct bound_box elaboration_box;
pippo = dim_vect;
obs = (struct Point *)G_calloc(pippo, sizeof(struct Point));
/* Reading points inside elaboration zone */
Vect_region_box(Elaboration, &elaboration_box);
npoints = 0;
nrows = Original->rows;
ncols = Original->cols;
if (Original->north > Elaboration->north)
startrow =
(Original->north - Elaboration->north) / Original->ns_res - 1;
else
startrow = 0;
if (Original->north > Elaboration->south) {
endrow = (Original->north - Elaboration->south) / Original->ns_res + 1;
if (endrow > nrows)
endrow = nrows;
}
else
endrow = nrows;
if (Elaboration->west > Original->west)
startcol = (Elaboration->west - Original->west) / Original->ew_res - 1;
else
startcol = 0;
if (Elaboration->east > Original->west) {
endcol = (Elaboration->east - Original->west) / Original->ew_res + 1;
if (endcol > ncols)
endcol = ncols;
}
else
endcol = ncols;
for (row = startrow; row < endrow; row++) {
for (col = startcol; col < endcol; col++) {
Segment_get(in_seg, &z, row, col);
if (!Rast_is_d_null_value(&z)) {
x = Rast_col_to_easting((double)(col) + 0.5, Original);
y = Rast_row_to_northing((double)(row) + 0.5, Original);
if (Vect_point_in_box(x, y, 0, &elaboration_box)) {
npoints++;
if (npoints >= pippo) {
pippo += dim_vect;
obs = (struct Point *)G_realloc(
(void *)obs,
(signed int)pippo * sizeof(struct Point));
}
/* Storing observation vector */
obs[npoints - 1].coordX = x;
obs[npoints - 1].coordY = y;
obs[npoints - 1].coordZ = z;
}
}
}
}
*num_points = npoints;
return obs;
}
/*------------------------------------------------------------------------------------------------*/
int P_Create_Aux2_Table(dbDriver *driver, char *tab_name)
{
dbTable *auxiliar_tab;
dbColumn *column;
auxiliar_tab = db_alloc_table(2);
db_set_table_name(auxiliar_tab, tab_name);
db_set_table_description(auxiliar_tab, "Intermediate interpolated values");
column = db_get_table_column(auxiliar_tab, 0);
db_set_column_name(column, "ID");
db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);
column = db_get_table_column(auxiliar_tab, 1);
db_set_column_name(column, "Interp");
db_set_column_sqltype(column, DB_SQL_TYPE_REAL);
if (db_create_table(driver, auxiliar_tab) == DB_OK) {
G_debug(1, _("<%s> created in database."), tab_name);
return TRUE;
}
else
G_warning(_("<%s> has not been created in database."), tab_name);
return FALSE;
}
/*------------------------------------------------------------------------------------------------*/
int P_Create_Aux4_Table(dbDriver *driver, char *tab_name)
{
dbTable *auxiliar_tab;
dbColumn *column;
auxiliar_tab = db_alloc_table(4);
db_set_table_name(auxiliar_tab, tab_name);
db_set_table_description(auxiliar_tab, "Intermediate interpolated values");
column = db_get_table_column(auxiliar_tab, 0);
db_set_column_name(column, "ID");
db_set_column_sqltype(column, DB_SQL_TYPE_INTEGER);
column = db_get_table_column(auxiliar_tab, 1);
db_set_column_name(column, "Interp");
db_set_column_sqltype(column, DB_SQL_TYPE_REAL);
column = db_get_table_column(auxiliar_tab, 2);
db_set_column_name(column, "X");
db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
column = db_get_table_column(auxiliar_tab, 3);
db_set_column_name(column, "Y");
db_set_column_sqltype(column, DB_SQL_TYPE_DOUBLE_PRECISION);
if (db_create_table(driver, auxiliar_tab) == DB_OK) {
G_debug(1, _("<%s> created in database."), tab_name);
return TRUE;
}
else
G_warning(_("<%s> has not been created in database."), tab_name);
return FALSE;
}
/*------------------------------------------------------------------------------------------------*/
int P_Drop_Aux_Table(dbDriver *driver, char *tab_name)
{
dbString drop;
db_init_string(&drop);
db_append_string(&drop, "drop table ");
db_append_string(&drop, tab_name);
return db_execute_immediate(driver, &drop);
}
/*---------------------------------------------------------------------------------------*/
void P_Aux_to_Raster(double **matrix, int fd)
{
int ncols, col, nrows, row;
void *ptr, *raster;
nrows = Rast_window_rows();
ncols = Rast_window_cols();
raster = Rast_allocate_buf(DCELL_TYPE);
for (row = 0; row < nrows; row++) {
G_percent(row, nrows, 2);
Rast_set_d_null_value(raster, ncols);
for (col = 0, ptr = raster; col < ncols;
col++, ptr = G_incr_void_ptr(ptr, Rast_cell_size(DCELL_TYPE))) {
Rast_set_d_value(ptr, (DCELL)(matrix[row][col]), DCELL_TYPE);
}
Rast_put_d_row(fd, raster);
}
G_percent(row, nrows, 2);
}
/*------------------------------------------------------------------------------------------------*/
void P_Aux_to_Vector(struct Map_info *Map G_UNUSED, struct Map_info *Out,
dbDriver *driver, char *tab_name)
{
int more, type;
double coordX, coordY, coordZ;
struct line_pnts *point;
struct line_cats *cat;
dbTable *table;
dbColumn *column;
dbValue *value;
dbCursor cursor;
dbString sql;
char buf[1024];
point = Vect_new_line_struct();
cat = Vect_new_cats_struct();
db_init_string(&sql);
db_zero_string(&sql);
snprintf(buf, sizeof(buf),
"select ID, X, Y, sum(Interp) from %s group by ID, X, Y",
tab_name);
db_append_string(&sql, buf);
db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL);
while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
table = db_get_cursor_table(&cursor);
column = db_get_table_column(table, 0);
type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
if (type == DB_C_TYPE_INT)
value = db_get_column_value(column);
else
continue;
column = db_get_table_column(table, 1);
type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
if (type == DB_C_TYPE_DOUBLE)
value = db_get_column_value(column);
else
continue;
coordZ = db_get_value_double(value);
column = db_get_table_column(table, 2);
type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
if (type == DB_C_TYPE_DOUBLE)
value = db_get_column_value(column);
else
continue;
coordX = db_get_value_double(value);
column = db_get_table_column(table, 3);
type = db_sqltype_to_Ctype(db_get_column_sqltype(column));
if (type == DB_C_TYPE_DOUBLE)
value = db_get_column_value(column);
else
continue;
coordY = db_get_value_double(value);
Vect_copy_xyz_to_pnts(point, &coordX, &coordY, &coordZ, 1);
Vect_reset_cats(cat);
Vect_cat_set(cat, 1, 1);
Vect_write_line(Out, GV_POINT, point, cat);
}
return;
}
/*! DEFINITION OF THE SUBZONES
5: inside Overlap region
all others: inside General region but outside Overlap region
---------------------------------
| | | | | | | |
---------------------------------
| | | | | | | |
| | | | | | | |
| | | | | | | |
---------------------------------
| | |4| 3 |3| | |
---------------------------------
| | | | | | | |
| | |2| 5 |1| | |
| | | | | | | |
---------------------------------
| | |2| 1 |1| | |
---------------------------------
| | | | | | | |
| | | | | | | |
| | | | | | | |
---------------------------------
| | | | | | | |
---------------------------------
*/
G_incr_void_ptr
void * G_incr_void_ptr(const void *ptr, size_t size)
Advance void pointer.
Definition alloc.c:184
NULL
#define NULL
Definition ccmath.h:32
TRUE
#define TRUE
Definition dbfopen.c:56
FALSE
#define FALSE
Definition dbfopen.c:55
G_debug
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
G_get_window
void G_get_window(struct Cell_head *window)
Get the current region.
Definition get_window.c:47
G_warning
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
G_get_set_window
void G_get_set_window(struct Cell_head *window)
Get the current working window (region).
Definition gis/set_window.c:35
GENERAL_ROW
#define GENERAL_ROW
Definition lidar.h:40
P_BILINEAR
#define P_BILINEAR
Definition lidar.h:64
LAST_COLUMN
#define LAST_COLUMN
Definition lidar.h:45
FIRST_ROW
#define FIRST_ROW
Definition lidar.h:42
FIRST_COLUMN
#define FIRST_COLUMN
Definition lidar.h:44
LAST_ROW
#define LAST_ROW
Definition lidar.h:43
CONTOUR
#define CONTOUR
Definition lidar.h:39
GENERAL_COLUMN
#define GENERAL_COLUMN
Definition lidar.h:41
P_BICUBIC
#define P_BICUBIC
Definition lidar.h:65
G_percent
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition percent.c:61
Segment_get
int Segment_get(SEGMENT *SEG, void *buf, off_t row, off_t col)
Get value from segment file.
Definition segment/get.c:36
Point
Definition lidar.h:77
Point::cat
int cat
Definition lidar.h:82
Point::lineID
int lineID
Definition lidar.h:81
Point::coordY
double coordY
Definition lidar.h:79
Point::coordX
double coordX
Definition lidar.h:78
Point::coordZ
double coordZ
Definition lidar.h:80
Reg_dimens
Definition lidar.h:69
Reg_dimens::edge_v
double edge_v
Definition lidar.h:71
Reg_dimens::edge_h
double edge_h
Definition lidar.h:70
Reg_dimens::ew_size
double ew_size
Definition lidar.h:74
Reg_dimens::overlap
double overlap
Definition lidar.h:72
Reg_dimens::sn_size
double sn_size
Definition lidar.h:73
driver
Definition driver.h:27
x
#define x
P_Aux_to_Vector
void P_Aux_to_Vector(struct Map_info *Map, struct Map_info *Out, dbDriver *driver, char *tab_name)
Definition zones.c:565
P_Aux_to_Raster
void P_Aux_to_Raster(double **matrix, int fd)
Definition zones.c:540
P_Mean_Calc
double P_Mean_Calc(struct Cell_head *Elaboration, struct Point *obs, int npoints)
Definition zones.c:231
P_estimate_splinestep
double P_estimate_splinestep(struct Map_info *Map, double *dens, double *dist)
Definition zones.c:259
P_Read_Vector_Region_Map
struct Point * P_Read_Vector_Region_Map(struct Map_info *Map, struct Cell_head *Elaboration, int *num_points, int dim_vect, int layer)
Definition zones.c:322
P_get_edge
int P_get_edge(int interpolator, struct Reg_dimens *dim, double pe, double pn)
Definition zones.c:196
P_zero_dim
void P_zero_dim(struct Reg_dimens *dim)
Definition zones.c:9
P_get_BandWidth
int P_get_BandWidth(int interpolator, int nsplines)
Definition zones.c:218
P_Read_Raster_Region_Map
struct Point * P_Read_Raster_Region_Map(SEGMENT *in_seg, struct Cell_head *Elaboration, struct Cell_head *Original, int *num_points, int dim_vect)
Definition zones.c:387
P_Create_Aux2_Table
int P_Create_Aux2_Table(dbDriver *driver, char *tab_name)
Definition zones.c:465
P_set_dim
int P_set_dim(struct Reg_dimens *dim, double pe, double pn, int *nsplx, int *nsply)
Definition zones.c:117
P_Drop_Aux_Table
int P_Drop_Aux_Table(dbDriver *driver, char *tab_name)
Definition zones.c:529
P_Create_Aux4_Table
int P_Create_Aux4_Table(dbDriver *driver, char *tab_name)
Definition zones.c:493
P_set_regions
int P_set_regions(struct Cell_head *Elaboration, struct bound_box *General, struct bound_box *Overlap, struct Reg_dimens dim, int type)
Definition zones.c:53