GRASS GIS 7 Programmer's Manual  7.5.svn(2017)-r71759
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
point2d.c
Go to the documentation of this file.
1 /*!
2  * \file point2d.c
3  *
4  * \author
5  * Lubos Mitas (original program and various modifications)
6  *
7  * \author
8  * H. Mitasova,
9  * I. Kosinovsky,
10  * D. Gerdes,
11  * D. McCauley
12  * (GRASS4.1 version of the program and GRASS4.2 modifications)
13  *
14  * \author modified by McCauley in August 1995
15  * \author modified by Mitasova in August 1995, Nov. 1996
16  *
17  * \copyright
18  * (C) 1993-2006 by Helena Mitasova and the GRASS Development Team
19  *
20  * \copyright
21  * This program is free software under the
22  * GNU General Public License (>=v2).
23  * Read the file COPYING that comes with GRASS for details.
24  */
25 
26 
27 #include <stdio.h>
28 #include <math.h>
29 #include <unistd.h>
30 #include <grass/gis.h>
31 #include <grass/vector.h>
32 #include <grass/dbmi.h>
33 
34 #define POINT2D_C
35 #include <grass/interpf.h>
36 
37 /* needed for AIX */
38 #ifdef hz
39 #undef hz
40 #endif
41 
42 /*!
43  * Checks if interpolating function interp() evaluates correct z-values at
44  * given points. If smoothing is used calculate the maximum error caused
45  * by smoothing.
46  *
47  * *ertot* is a RMS deviation of the interpolated surface.
48  *
49  * \todo
50  * Alternative description:
51  * ...calculate the maximum and RMS deviation caused by smoothing.
52  */
54  struct quaddata *data, /*!< current region */
55  double *b, /*!< solution of linear equations */
56  double *ertot, /*!< total error */
57  double zmin, /*!< min z-value */
58  double dnorm,
59  struct triple skip_point
60  )
61 {
62  int n_points = data->n_points; /* number of points */
63  struct triple *points = data->points; /* points for interpolation */
64  double east = data->xmax;
65  double west = data->x_orig;
66  double north = data->ymax;
67  double south = data->y_orig;
68  double rfsta2, errmax, h, xx, yy, r2, hz, zz, err, xmm, ymm, r;
69  double skip_err;
70  int n1, mm, m, cat;
71  double fstar2;
72  int inside;
73 
74  /* Site *site; */
75  char buf[1024];
76 
77 
78  /* if ((site = G_site_new_struct (-1, 2, 0, 1)) == NULL)
79  G_fatal_error ("Memory error for site struct"); */
80 
81  fstar2 = params->fi * params->fi / 4.;
82  errmax = .0;
83  n1 = n_points + 1;
84  for (mm = 1; mm <= n_points; mm++) {
85  h = b[0];
86  for (m = 1; m <= n_points; m++) {
87  xx = points[mm - 1].x - points[m - 1].x;
88  yy = points[mm - 1].y - points[m - 1].y;
89  r2 = yy * yy + xx * xx;
90  if (r2 != 0.) {
91  rfsta2 = fstar2 * r2;
92  r = r2;
93  h = h + b[m] * params->interp(r, params->fi);
94  }
95  }
96  /* modified by helena january 1997 - normalization of z was
97  removed from segm2d.c and interp2d.c
98  hz = (h * dnorm) + zmin;
99  zz = (points[mm - 1].z * dnorm) + zmin;
100  */
101  hz = h + zmin;
102  zz = points[mm - 1].z + zmin;
103  err = hz - zz;
104  xmm = points[mm - 1].x * dnorm + params->x_orig + west;
105  ymm = points[mm - 1].y * dnorm + params->y_orig + south;
106  if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig)
107  && (ymm >= south + params->y_orig) &&
108  (ymm <= north + params->y_orig))
109  inside = 1;
110  else
111  inside = 0;
112 
113  if (params->fddevi != NULL) {
114 
115  if (inside) { /* if the point is inside the region */
118 
119  Vect_append_point(Pnts, xmm, ymm, zz);
120  cat = count;
121  Vect_cat_set(Cats2, 1, cat);
123 
125  sprintf(buf, "insert into %s values ( %d ", ff->table, cat);
126  db_append_string(&sql2, buf);
127 
128  sprintf(buf, ", %f", err);
129  db_append_string(&sql2, buf);
130  db_append_string(&sql2, ")");
131  G_debug(3, "IL_check_at_points_2d: %s", db_get_string(&sql2));
132 
136  G_fatal_error("Cannot insert new row: %s",
137  db_get_string(&sql2));
138  }
139  count++;
140 
141  }
142  }
143  (*ertot) += err * err;
144  }
145 
146  /* cv stuff */
147  if (params->cv) {
148  h = b[0];
149  for (m = 1; m <= n_points - 1; m++) {
150  xx = points[m - 1].x - skip_point.x;
151  yy = points[m - 1].y - skip_point.y;
152  r2 = yy * yy + xx * xx;
153  if (r2 != 0.) {
154  rfsta2 = fstar2 * r2;
155  r = r2;
156  h = h + b[m] * params->interp(r, params->fi);
157  }
158  }
159  hz = h + zmin;
160  zz = skip_point.z + zmin;
161  skip_err = hz - zz;
162  xmm = skip_point.x * dnorm + params->x_orig + west;
163  ymm = skip_point.y * dnorm + params->y_orig + south;
164 
165  if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig)
166  && (ymm >= south + params->y_orig) &&
167  (ymm <= north + params->y_orig))
168  inside = 1;
169  else
170  inside = 0;
171 
172  if (inside) { /* if the point is inside the region */
175 
176  Vect_append_point(Pnts, xmm, ymm, zz);
177  cat = count;
178  Vect_cat_set(Cats2, 1, cat);
180 
182  sprintf(buf, "insert into %s values ( %d ", ff->table, cat);
183  db_append_string(&sql2, buf);
184 
185  sprintf(buf, ", %f", skip_err);
186  db_append_string(&sql2, buf);
187  db_append_string(&sql2, ")");
188  G_debug(3, "IL_check_at_points_2d: %s", db_get_string(&sql2));
189 
193  G_fatal_error("Cannot insert new row: %s",
194  db_get_string(&sql2));
195  }
196  count++;
197  }
198  } /* cv */
199 
200 
201  return 1;
202 }
void db_zero_string(dbString *x)
Zero string.
Definition: string.c:79
interp_fn * interp
Definition: interpf.h:101
CELL cat
Definition: raster3d/cats.c:82
double y_orig
Definition: dataquad.h:51
dbString sql2
off_t Vect_write_line(struct Map_info *Map, int type, const struct line_pnts *points, const struct line_cats *cats)
Writes a new feature.
int db_shutdown_driver(dbDriver *driver)
Closedown the driver, and free the driver structure.
Definition: shutdown.c:36
int db_close_database(dbDriver *driver)
Close database connection.
Definition: c_closedb.c:26
char * table
Name of DB table.
Definition: dig_structs.h:155
double z
Definition: dataquad.h:44
int count
double y_orig
Definition: interpf.h:87
int Vect_reset_cats(struct line_cats *Cats)
Reset category structure to make sure cats structure is clean to be re-used.
#define NULL
Definition: ccmath.h:32
double x_orig
Definition: dataquad.h:50
int IL_check_at_points_2d(struct interp_params *, struct quaddata *, double *, double *, double, double, struct triple)
Definition: point2d.c:53
#define GV_POINT
Feature types used in memory on run time (may change)
Definition: dig_defines.h:182
char * db_get_string(const dbString *x)
Get string.
Definition: string.c:140
int db_execute_immediate(dbDriver *driver, dbString *SQLstatement)
Execute SQL statements.
Definition: c_execute.c:27
int Vect_append_point(struct line_pnts *Points, double x, double y, double z)
Appends one point to the end of a line.
Definition: line.c:149
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:159
struct triple * points
Definition: dataquad.h:57
int db_append_string(dbString *x, const char *s)
Append string to dbString.
Definition: string.c:205
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:220
int Vect_cat_set(struct line_cats *Cats, int field, int cat)
Add new field/cat to category structure if doesn&#39;t exist yet.
FILE * fddevi
Definition: interpf.h:95
dbDriver * driver2
double b
Definition: r_raster.c:39
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
struct field_info * ff
double x
Definition: dataquad.h:42
double fi
Definition: interpf.h:80
double x_orig
Definition: interpf.h:87
double ymax
Definition: dataquad.h:53
void Vect_reset_line(struct line_pnts *Points)
Reset line.
Definition: line.c:130
double xmax
Definition: dataquad.h:52
struct line_pnts * Pnts
double y
Definition: dataquad.h:43
struct Map_info Map2
double r
Definition: r_raster.c:39
struct line_cats * Cats2
#define DB_OK
Definition: dbmi.h:71
int n_points
Definition: dataquad.h:56