Loading...
Searching...
No Matches
vinput2d.c
Go to the documentation of this file.
/*!
* \file vinput2d.c
*
* \author
* Written by H. Mitasova, I. Kosinovsky, D. Gerdes Fall 1993
* University of Illinois
* US Army Construction Engineering Research Lab
*
* \author
* Mitasova (University of Illinois),
* I. Kosinovsky, (USA-CERL), and D.Gerdes (USA-CERL)
*
* \author modified by McCauley in August 1995
* \author modified by Mitasova in August 1995
* \author modofied by Mitasova in Nov 1999 (dmax fix)
*
* \copyright
* (C) 1993-1999 by Helena Mitasova and the GRASS Development Team
*
* \copyright
* This program is free software under the
* GNU General Public License (>=v2).
* Read the file COPYING that comes with GRASS
* for details.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <grass/bitmap.h>
#include <grass/linkm.h>
#include <grass/gis.h>
#include <grass/dbmi.h>
#include <grass/vector.h>
#include <grass/glocale.h>
#include <grass/interpf.h>
/*!
* Insert into a quad tree
*
* Inserts input data inside the region into a quad tree. Also translates
* data. Returns number of segments in the quad tree.
*
* As z values may be used (in *Map*):
* - z coordinates in 3D file -> field = 0
* - categories -> field > 0, zcol = NULL
* - attributes -> field > 0, zcol != NULL
*/
int IL_vector_input_data_2d(
struct interp_params *params, /*!< interpolation parameters */
struct Map_info *Map, /*!< input vector map */
int field, /*!< category field number */
char *zcol, /*!< name of the column containing z values */
char *scol, /*!< name of the column containing smooth values */
struct tree_info *info, /*!< quadtree info */
double *xmin, double *xmax, double *ymin, double *ymax, double *zmin,
double *zmax, int *n_points, /*!< number of points used for interpolation */
double *dmax /*!< max distance between points */
)
{
double dmax2; /* max distance between points squared */
double c1, c2, c3, c4;
int i, k = 0;
double ns_res, ew_res;
int npoint, OUTRANGE;
int totsegm;
struct quaddata *data = (struct quaddata *)info->root->data;
double xprev, yprev, zprev, x1, y1, z1, d1, xt, yt, z, sm;
struct line_pnts *Points;
struct line_cats *Cats;
int times, j1, ltype, cat, zctype = 0, sctype = 0;
struct field_info *Fi;
dbDriver *driver;
dbHandle handle;
dbString stmt;
dbCatValArray zarray, sarray;
OUTRANGE = 0;
npoint = 0;
G_debug(2, "IL_vector_input_data_2d(): field = %d, zcol = %s, scol = %s",
field, zcol, scol);
ns_res = (data->ymax - data->y_orig) / data->n_rows;
ew_res = (data->xmax - data->x_orig) / data->n_cols;
dmax2 = *dmax * *dmax;
Points = Vect_new_line_struct(); /* init line_pnts struct */
Cats = Vect_new_cats_struct();
if (field == 0 && !Vect_is_3d(Map))
G_fatal_error(_("Vector map <%s> is not 3D"), Vect_get_full_name(Map));
if (field > 0 && zcol != NULL) { /* open db driver */
G_verbose_message(_("Loading data from attribute table ..."));
Fi = Vect_get_field(Map, field);
if (Fi == NULL)
G_fatal_error(_("Database connection not defined for layer %d"),
field);
G_debug(3, " driver = %s database = %s table = %s", Fi->driver,
Fi->database, Fi->table);
db_init_handle(&handle);
db_init_string(&stmt);
driver = db_start_driver(Fi->driver);
db_set_handle(&handle, Fi->database, NULL);
if (db_open_database(driver, &handle) != DB_OK)
G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
Fi->database, Fi->driver);
zctype = db_column_Ctype(driver, Fi->table, zcol);
G_debug(3, " zcol C type = %d", zctype);
if (zctype == -1)
G_fatal_error(_("Column <%s> not found"), zcol);
if (zctype != DB_C_TYPE_INT && zctype != DB_C_TYPE_DOUBLE)
G_fatal_error(_("Data type of column <%s> must be numeric"), zcol);
db_CatValArray_init(&zarray);
G_debug(3, "RST SQL WHERE: %s", params->wheresql);
db_select_CatValArray(driver, Fi->table, Fi->key, zcol,
params->wheresql, &zarray);
if (scol != NULL) {
sctype = db_column_Ctype(driver, Fi->table, scol);
G_debug(3, " scol C type = %d", sctype);
if (sctype == -1)
G_fatal_error(_("Column <%s> not found"), scol);
if (sctype != DB_C_TYPE_INT && sctype != DB_C_TYPE_DOUBLE)
G_fatal_error(_("Data type of column <%s> must be numeric"),
scol);
db_CatValArray_init(&sarray);
db_select_CatValArray(driver, Fi->table, Fi->key, scol,
params->wheresql, &sarray);
}
db_close_database_shutdown_driver(driver);
}
/* Lines without nodes */
G_message(_("Reading features from vector map ..."));
sm = 0;
while ((ltype = Vect_read_next_line(Map, Points, Cats)) != -2) {
if (!(ltype & (GV_POINT | GV_LINE | GV_BOUNDARY)))
continue;
if (field > 0) { /* use cat or attribute */
Vect_cat_get(Cats, field, &cat);
if (zcol == NULL) { /* use categories */
z = (double)cat;
}
else { /* read att from db */
int ret, intval;
if (zctype == DB_C_TYPE_INT) {
ret = db_CatValArray_get_value_int(&zarray, cat, &intval);
z = intval;
}
else { /* DB_C_TYPE_DOUBLE */
ret = db_CatValArray_get_value_double(&zarray, cat, &z);
}
if (ret != DB_OK) {
if (params->wheresql != NULL)
/* G_message(_("Database record for cat %d not used due
* to SQL statement")); */
/* do nothing in this case to not confuse user. Or
* implement second cat list */
;
else
G_warning(_("Database record for cat %d not found"),
cat);
continue;
}
if (scol != NULL) {
if (sctype == DB_C_TYPE_INT) {
ret =
db_CatValArray_get_value_int(&sarray, cat, &intval);
sm = intval;
}
else { /* DB_C_TYPE_DOUBLE */
ret =
db_CatValArray_get_value_double(&sarray, cat, &sm);
}
if (sm < 0.0)
G_fatal_error(_("Negative value of smoothing detected: "
"sm must be >= 0"));
}
G_debug(5, " z = %f sm = %f", z, sm);
}
}
/* Insert all points including nodes (end points) */
for (i = 0; i < Points->n_points; i++) {
if (field == 0)
z = Points->z[i];
process_point(Points->x[i], Points->y[i], z, sm, info,
params->zmult, xmin, xmax, ymin, ymax, zmin, zmax,
&npoint, &OUTRANGE, &k);
}
/* Check all segments */
xprev = Points->x[0];
yprev = Points->y[0];
zprev = Points->z[0];
for (i = 1; i < Points->n_points; i++) {
/* compare the distance between current and previous */
x1 = Points->x[i];
y1 = Points->y[i];
z1 = Points->z[i];
xt = x1 - xprev;
yt = y1 - yprev;
d1 = (xt * xt + yt * yt);
if ((d1 > dmax2) && (dmax2 != 0.)) {
times = (int)(d1 / dmax2 + 0.5);
for (j1 = 0; j1 < times; j1++) {
xt = x1 - j1 * ((x1 - xprev) / times);
yt = y1 - j1 * ((y1 - yprev) / times);
if (field == 0)
z = z1 - j1 * ((z1 - zprev) / times);
process_point(xt, yt, z, sm, info, params->zmult, xmin,
xmax, ymin, ymax, zmin, zmax, &npoint,
&OUTRANGE, &k);
}
}
xprev = x1;
yprev = y1;
zprev = z1;
}
}
if (field > 0 && zcol != NULL)
db_CatValArray_free(&zarray);
if (scol != NULL) {
db_CatValArray_free(&sarray);
}
c1 = *xmin - data->x_orig;
c2 = data->xmax - *xmax;
c3 = *ymin - data->y_orig;
c4 = data->ymax - *ymax;
if ((c1 > 5 * ew_res) || (c2 > 5 * ew_res) || (c3 > 5 * ns_res) ||
(c4 > 5 * ns_res)) {
static int once = 0;
if (!once) {
once = 1;
G_warning(_("Strip exists with insufficient data"));
}
}
totsegm = translate_quad(info->root, data->x_orig, data->y_orig, *zmin, 4);
if (!totsegm)
return 0;
data->x_orig = 0;
data->y_orig = 0;
/* G_read_vector_timestamp(name,mapset,ts); */
if (OUTRANGE > 0)
G_warning(_("There are points outside specified 2D/3D region - %d "
"points ignored"),
OUTRANGE);
if (npoint > 0)
G_important_message(_("Ignoring %d points (too dense)"), npoint);
npoint = k - npoint - OUTRANGE;
if (npoint < params->kmin) {
if (npoint != 0) {
G_warning(_("%d points given for interpolation (after thinning) is "
"less than given NPMIN=%d"),
npoint, params->kmin);
params->kmin = npoint;
}
else {
G_warning(_("Zero points in the given region"));
return -1;
}
}
if (npoint > params->KMAX2 && params->kmin <= params->kmax) {
G_warning(
_("Segmentation parameters set to invalid values: npmin= %d, "
"segmax= %d "
"for smooth connection of segments, npmin > segmax (see manual)"),
params->kmin, params->kmax);
return -1;
}
if (npoint < params->KMAX2 && params->kmax != params->KMAX2)
G_warning(_("There are less than %d points for interpolation. No "
"segmentation is necessary, to run the program faster set "
"segmax=%d (see manual)"),
params->KMAX2, params->KMAX2);
G_verbose_message(_("Number of points from vector map %d"), k);
G_verbose_message(_("Number of points outside of 2D/3D region %d"),
OUTRANGE);
G_verbose_message(_("Number of points being used %d"), npoint);
*n_points = npoint;
return (totsegm);
}
int process_point(double x, double y, double z, double sm,
struct tree_info *info, /* quadtree info */
double zmult, /* multiplier for z-values */
double *xmin, double *xmax, double *ymin, double *ymax,
double *zmin, double *zmax, int *npoint, int *OUTRANGE,
int *total)
{
struct triple *point;
double c1, c2, c3, c4;
int a;
static int first_time = 1;
struct quaddata *data = (struct quaddata *)info->root->data;
(*total)++;
z = z * zmult;
c1 = x - data->x_orig;
c2 = data->xmax - x;
c3 = y - data->y_orig;
c4 = data->ymax - y;
if (!((c1 >= 0) && (c2 >= 0) && (c3 >= 0) && (c4 >= 0))) {
if (!(*OUTRANGE)) {
G_warning(_("Some points outside of region (ignored)"));
}
(*OUTRANGE)++;
}
else {
if (!(point = quad_point_new(x, y, z, sm))) {
G_warning(_("Unable to allocate memory"));
return -1;
}
a = MT_insert(point, info, info->root, 4);
if (a == 0) {
(*npoint)++;
}
if (a < 0) {
G_warning(_("Unable to insert %f,%f,%f a = %d"), x, y, z, a);
return -1;
}
free(point);
if (first_time) {
first_time = 0;
*xmin = x;
*ymin = y;
*zmin = z;
*xmax = x;
*ymax = y;
*zmax = z;
}
*xmin = amin1(*xmin, x);
*ymin = amin1(*ymin, y);
*zmin = amin1(*zmin, z);
*xmax = amax1(*xmax, x);
*ymax = amax1(*ymax, y);
*zmax = amax1(*zmax, z);
}
return 1;
}
NULL
#define NULL
Definition ccmath.h:32
quad_point_new
struct triple * quad_point_new(double x, double y, double z, double sm)
Definition dataquad.c:36
G_debug
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
G_verbose_message
void G_verbose_message(const char *msg,...)
Print a message to stderr but only if module is in verbose mode.
Definition gis/error.c:108
G_important_message
void G_important_message(const char *msg,...)
Print a message to stderr even in brief mode (verbosity=1).
Definition gis/error.c:130
G_fatal_error
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
G_message
void G_message(const char *msg,...)
Print a message to stderr.
Definition gis/error.c:89
G_warning
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
translate_quad
int translate_quad(struct multtree *tree, double numberx, double numbery, double numberz, int n_leafs)
Definition input2d.c:90
amin1
double amin1(double, double)
Definition minmax.c:65
amax1
double amax1(double, double)
Definition minmax.c:52
MT_insert
int MT_insert(struct triple *point, struct tree_info *info, struct multtree *tree, int n_leafs)
Definition qtree.c:103
driver
Definition driver.h:27
interp_params
Definition interpf.h:70
interp_params::zmult
double zmult
Definition interpf.h:72
interp_params::wheresql
const char * wheresql
Definition interpf.h:142
interp_params::kmin
int kmin
Definition interpf.h:84
interp_params::kmax
int kmax
Definition interpf.h:87
interp_params::KMAX2
int KMAX2
Definition interpf.h:99
multtree::data
struct quaddata * data
Definition qtree.h:53
quaddata
Definition dataquad.h:45
quaddata::ymax
double ymax
Definition dataquad.h:49
quaddata::xmax
double xmax
Definition dataquad.h:48
tree_info
Definition qtree.h:45
tree_info::root
struct multtree * root
Definition qtree.h:49
triple
Definition dataquad.h:38
process_point
int process_point(double x, double y, double z, double sm, struct tree_info *info, double zmult, double *xmin, double *xmax, double *ymin, double *ymax, double *zmin, double *zmax, int *npoint, int *OUTRANGE, int *total)
Definition vinput2d.c:306
IL_vector_input_data_2d
int IL_vector_input_data_2d(struct interp_params *params, struct Map_info *Map, int field, char *zcol, char *scol, struct tree_info *info, double *xmin, double *xmax, double *ymin, double *ymax, double *zmin, double *zmax, int *n_points, double *dmax)
Definition vinput2d.c:50
x
#define x