Loading...
Searching...
No Matches
point2d.c
Go to the documentation of this file.
/*!
* \file point2d.c
*
* \author
* Lubos Mitas (original program and various modifications)
*
* \author
* H. Mitasova,
* I. Kosinovsky,
* D. Gerdes,
* D. McCauley
* (GRASS4.1 version of the program and GRASS4.2 modifications)
*
* \author modified by McCauley in August 1995
* \author modified by Mitasova in August 1995, Nov. 1996
*
* \copyright
* (C) 1993-2006 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 <math.h>
#include <unistd.h>
#include <grass/gis.h>
#include <grass/vector.h>
#include <grass/dbmi.h>
#define POINT2D_C
#include <grass/interpf.h>
/* needed for AIX */
#ifdef hz
#undef hz
#endif
/*!
* Checks if interpolating function interp() evaluates correct z-values at
* given points. If smoothing is used calculate the maximum error caused
* by smoothing.
*
* *ertot* is a RMS deviation of the interpolated surface.
*
* \todo
* Alternative description:
* ...calculate the maximum and RMS deviation caused by smoothing.
*/
int IL_check_at_points_2d(struct interp_params *params,
struct quaddata *data, /*!< current region */
double *b, /*!< solution of linear equations */
double *ertot, /*!< total error */
double zmin, /*!< min z-value */
double dnorm, struct triple *skip_point)
{
int n_points = data->n_points; /* number of points */
struct triple *points = data->points; /* points for interpolation */
struct triple point_writeout;
double east = data->xmax;
double west = data->x_orig;
double north = data->ymax;
double south = data->y_orig;
double /* rfsta2, errmax, */ h, xx, yy, r2, hz, zz, err, xmm, ymm, r;
double skip_err;
int /* n1, */ mm, m;
/* double fstar2; */
int inside;
/* if ((site = G_site_new_struct (-1, 2, 0, 1)) == NULL)
G_fatal_error ("Memory error for site struct"); */
/* fstar2 = params->fi * params->fi / 4.; */
/* errmax = .0; */
/* n1 = n_points + 1; */
for (mm = 1; mm <= n_points; mm++) {
h = b[0];
for (m = 1; m <= n_points; m++) {
xx = points[mm - 1].x - points[m - 1].x;
yy = points[mm - 1].y - points[m - 1].y;
r2 = yy * yy + xx * xx;
if (r2 != 0.) {
/* rfsta2 = fstar2 * r2; */
r = r2;
h = h + b[m] * params->interp(r, params->fi);
}
}
/* modified by helena january 1997 - normalization of z was
removed from segm2d.c and interp2d.c
hz = (h * dnorm) + zmin;
zz = (points[mm - 1].z * dnorm) + zmin;
*/
hz = h + zmin;
zz = points[mm - 1].z + zmin;
err = hz - zz;
xmm = points[mm - 1].x * dnorm + params->x_orig + west;
ymm = points[mm - 1].y * dnorm + params->y_orig + south;
if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig) &&
(ymm >= south + params->y_orig) && (ymm <= north + params->y_orig))
inside = 1;
else
inside = 0;
if (params->create_devi) {
if (inside) { /* if the point is inside the region */
point_writeout.x = xmm;
point_writeout.y = ymm;
point_writeout.z = zz;
IL_write_point_2d(point_writeout, err);
}
}
(*ertot) += err * err;
}
/* cv stuff */
if (params->cv) {
h = b[0];
for (m = 1; m <= n_points - 1; m++) {
xx = points[m - 1].x - skip_point->x;
yy = points[m - 1].y - skip_point->y;
r2 = yy * yy + xx * xx;
if (r2 != 0.) {
/* rfsta2 = fstar2 * r2; */
r = r2;
h = h + b[m] * params->interp(r, params->fi);
}
}
hz = h + zmin;
zz = skip_point->z + zmin;
skip_err = hz - zz;
xmm = skip_point->x * dnorm + params->x_orig + west;
ymm = skip_point->y * dnorm + params->y_orig + south;
if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig) &&
(ymm >= south + params->y_orig) && (ymm <= north + params->y_orig))
inside = 1;
else
inside = 0;
if (inside) { /* if the point is inside the region */
point_writeout.x = xmm;
point_writeout.y = ymm;
point_writeout.z = zz;
IL_write_point_2d(point_writeout, skip_err);
}
} /* cv */
return 1;
}
/*!
* \brief A function to write out point and deviation at point to database.
*
* \param point point to write out
* \param error deviation at point
*
* \return 1
*/
int IL_write_point_2d(struct triple point, double err)
{
char buf[1024];
Vect_reset_line(Pnts);
Vect_reset_cats(Cats2);
Vect_append_point(Pnts, point.x, point.y, point.z);
Vect_cat_set(Cats2, 1, count);
Vect_write_line(&Map2, GV_POINT, Pnts, Cats2);
db_zero_string(&sql2);
snprintf(buf, sizeof(buf), "insert into %s values ( %d ", ff->table, count);
db_append_string(&sql2, buf);
snprintf(buf, sizeof(buf), ", %f", err);
db_append_string(&sql2, buf);
db_append_string(&sql2, ")");
G_debug(3, "IL_check_at_points_2d: %s", db_get_string(&sql2));
if (db_execute_immediate(driver2, &sql2) != DB_OK) {
db_close_database(driver2);
db_shutdown_driver(driver2);
G_fatal_error("Cannot insert new row: %s", db_get_string(&sql2));
}
count++;
return 1;
}
G_debug
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
b
double b
Definition driver/set_window.c:5
r
double r
Definition driver/set_window.c:5
G_fatal_error
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
Map2
struct Map_info Map2
ff
struct field_info * ff
sql2
dbString sql2
driver2
dbDriver * driver2
Cats2
struct line_cats * Cats2
count
int count
Pnts
struct line_pnts * Pnts
IL_check_at_points_2d
int IL_check_at_points_2d(struct interp_params *params, struct quaddata *data, double *b, double *ertot, double zmin, double dnorm, struct triple *skip_point)
Definition point2d.c:52
IL_write_point_2d
int IL_write_point_2d(struct triple point, double err)
A function to write out point and deviation at point to database.
Definition point2d.c:163
interp_params
Definition interpf.h:70
interp_params::interp
interp_fn * interp
Definition interpf.h:136
interp_params::fi
double fi
Definition interpf.h:97
interp_params::x_orig
double x_orig
Definition interpf.h:112
interp_params::cv
int cv
Definition interpf.h:114
interp_params::y_orig
double y_orig
Definition interpf.h:112
interp_params::create_devi
bool create_devi
Definition interpf.h:126
quaddata
Definition dataquad.h:45
triple
Definition dataquad.h:38
triple::z
double z
Definition dataquad.h:41
triple::x
double x
Definition dataquad.h:39
triple::y
double y
Definition dataquad.h:40
err
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition symbol/read.c:216