GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-b2e2571e19
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 #include <stdio.h>
27 #include <math.h>
28 #include <unistd.h>
29 #include <grass/gis.h>
30 #include <grass/vector.h>
31 #include <grass/dbmi.h>
32 
33 #define POINT2D_C
34 #include <grass/interpf.h>
35 
36 /* needed for AIX */
37 #ifdef hz
38 #undef hz
39 #endif
40 
41 /*!
42  * Checks if interpolating function interp() evaluates correct z-values at
43  * given points. If smoothing is used calculate the maximum error caused
44  * by smoothing.
45  *
46  * *ertot* is a RMS deviation of the interpolated surface.
47  *
48  * \todo
49  * Alternative description:
50  * ...calculate the maximum and RMS deviation caused by smoothing.
51  */
53  struct quaddata *data, /*!< current region */
54  double *b, /*!< solution of linear equations */
55  double *ertot, /*!< total error */
56  double zmin, /*!< min z-value */
57  double dnorm, struct triple *skip_point)
58 {
59  int n_points = data->n_points; /* number of points */
60  struct triple *points = data->points; /* points for interpolation */
61  struct triple point_writeout;
62  double east = data->xmax;
63  double west = data->x_orig;
64  double north = data->ymax;
65  double south = data->y_orig;
66  double /* rfsta2, errmax, */ h, xx, yy, r2, hz, zz, err, xmm, ymm, r;
67  double skip_err;
68  int /* n1, */ mm, m;
69 
70  /* double fstar2; */
71  int inside;
72 
73  /* if ((site = G_site_new_struct (-1, 2, 0, 1)) == NULL)
74  G_fatal_error ("Memory error for site struct"); */
75 
76  /* fstar2 = params->fi * params->fi / 4.; */
77  /* errmax = .0; */
78  /* n1 = n_points + 1; */
79  for (mm = 1; mm <= n_points; mm++) {
80  h = b[0];
81  for (m = 1; m <= n_points; m++) {
82  xx = points[mm - 1].x - points[m - 1].x;
83  yy = points[mm - 1].y - points[m - 1].y;
84  r2 = yy * yy + xx * xx;
85  if (r2 != 0.) {
86  /* rfsta2 = fstar2 * r2; */
87  r = r2;
88  h = h + b[m] * params->interp(r, params->fi);
89  }
90  }
91  /* modified by helena january 1997 - normalization of z was
92  removed from segm2d.c and interp2d.c
93  hz = (h * dnorm) + zmin;
94  zz = (points[mm - 1].z * dnorm) + zmin;
95  */
96  hz = h + zmin;
97  zz = points[mm - 1].z + zmin;
98  err = hz - zz;
99  xmm = points[mm - 1].x * dnorm + params->x_orig + west;
100  ymm = points[mm - 1].y * dnorm + params->y_orig + south;
101  if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig) &&
102  (ymm >= south + params->y_orig) && (ymm <= north + params->y_orig))
103  inside = 1;
104  else
105  inside = 0;
106 
107  if (params->create_devi) {
108 
109  if (inside) { /* if the point is inside the region */
110  point_writeout.x = xmm;
111  point_writeout.y = ymm;
112  point_writeout.z = zz;
113  IL_write_point_2d(point_writeout, err);
114  }
115  }
116  (*ertot) += err * err;
117  }
118 
119  /* cv stuff */
120  if (params->cv) {
121  h = b[0];
122  for (m = 1; m <= n_points - 1; m++) {
123  xx = points[m - 1].x - skip_point->x;
124  yy = points[m - 1].y - skip_point->y;
125  r2 = yy * yy + xx * xx;
126  if (r2 != 0.) {
127  /* rfsta2 = fstar2 * r2; */
128  r = r2;
129  h = h + b[m] * params->interp(r, params->fi);
130  }
131  }
132  hz = h + zmin;
133  zz = skip_point->z + zmin;
134  skip_err = hz - zz;
135  xmm = skip_point->x * dnorm + params->x_orig + west;
136  ymm = skip_point->y * dnorm + params->y_orig + south;
137 
138  if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig) &&
139  (ymm >= south + params->y_orig) && (ymm <= north + params->y_orig))
140  inside = 1;
141  else
142  inside = 0;
143 
144  if (inside) { /* if the point is inside the region */
145  point_writeout.x = xmm;
146  point_writeout.y = ymm;
147  point_writeout.z = zz;
148  IL_write_point_2d(point_writeout, skip_err);
149  }
150  } /* cv */
151 
152  return 1;
153 }
154 
155 /*!
156  * \brief A function to write out point and deviation at point to database.
157  *
158  * \param point point to write out
159  * \param error deviation at point
160  *
161  * \return 1
162  */
163 
164 int IL_write_point_2d(struct triple point, double err)
165 {
166 
167  char buf[1024];
168 
171 
172  Vect_append_point(Pnts, point.x, point.y, point.z);
173  Vect_cat_set(Cats2, 1, count);
175 
177  sprintf(buf, "insert into %s values ( %d ", ff->table, count);
178  db_append_string(&sql2, buf);
179 
180  sprintf(buf, ", %f", err);
181  db_append_string(&sql2, buf);
182  db_append_string(&sql2, ")");
183  G_debug(3, "IL_check_at_points_2d: %s", db_get_string(&sql2));
184 
188  G_fatal_error("Cannot insert new row: %s", db_get_string(&sql2));
189  }
190  count++;
191 
192  return 1;
193 }
#define DB_OK
Definition: dbmi.h:71
int db_shutdown_driver(dbDriver *)
Closedown the driver, and free the driver structure.
Definition: shutdown.c:36
void db_zero_string(dbString *)
Zero string.
Definition: string.c:79
int db_execute_immediate(dbDriver *, dbString *)
Execute SQL statements.
Definition: c_execute.c:27
int db_close_database(dbDriver *)
Close database connection.
Definition: c_closedb.c:26
int db_append_string(dbString *, const char *)
Append string to dbString.
Definition: string.c:205
char * db_get_string(const dbString *)
Get string.
Definition: string.c:140
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
int Vect_reset_cats(struct line_cats *)
Reset category structure to make sure cats structure is clean to be re-used.
int Vect_cat_set(struct line_cats *, int, int)
Add new field/cat to category structure if doesn't exist yet.
off_t Vect_write_line(struct Map_info *, int, const struct line_pnts *, const struct line_cats *)
Writes a new feature.
void Vect_reset_line(struct line_pnts *)
Reset line.
Definition: line.c:129
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition: line.c:148
#define GV_POINT
Feature types used in memory on run time (may change)
Definition: dig_defines.h:183
struct Map_info Map2
struct field_info * ff
dbString sql2
dbDriver * driver2
struct line_cats * Cats2
int count
struct line_pnts * Pnts
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
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:164
double b
Definition: r_raster.c:39
double r
Definition: r_raster.c:39
char * table
Name of DB table.
Definition: dig_structs.h:151
interp_fn * interp
Definition: interpf.h:131
double fi
Definition: interpf.h:92
double x_orig
Definition: interpf.h:107
double y_orig
Definition: interpf.h:107
bool create_devi
Definition: interpf.h:121
double ymax
Definition: dataquad.h:49
double y_orig
Definition: dataquad.h:47
double x_orig
Definition: dataquad.h:46
struct triple * points
Definition: dataquad.h:53
int n_points
Definition: dataquad.h:52
double xmax
Definition: dataquad.h:48
double z
Definition: dataquad.h:41
double x
Definition: dataquad.h:39
double y
Definition: dataquad.h:40
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:216