GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
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 double east = data->xmax;
62 double west = data->x_orig;
63 double north = data->ymax;
64 double south = data->y_orig;
65 double /* rfsta2, errmax, */ h, xx, yy, r2, hz, zz, err, xmm, ymm, r;
66 double skip_err;
67 int /* n1, */ mm, m, cat;
68
69 /* double fstar2; */
70 int inside;
71
72 /* Site *site; */
73 char buf[1024];
74
75 /* if ((site = G_site_new_struct (-1, 2, 0, 1)) == NULL)
76 G_fatal_error ("Memory error for site struct"); */
77
78 /* fstar2 = params->fi * params->fi / 4.; */
79 /* errmax = .0; */
80 /* n1 = n_points + 1; */
81 for (mm = 1; mm <= n_points; mm++) {
82 h = b[0];
83 for (m = 1; m <= n_points; m++) {
84 xx = points[mm - 1].x - points[m - 1].x;
85 yy = points[mm - 1].y - points[m - 1].y;
86 r2 = yy * yy + xx * xx;
87 if (r2 != 0.) {
88 /* rfsta2 = fstar2 * r2; */
89 r = r2;
90 h = h + b[m] * params->interp(r, params->fi);
91 }
92 }
93 /* modified by helena january 1997 - normalization of z was
94 removed from segm2d.c and interp2d.c
95 hz = (h * dnorm) + zmin;
96 zz = (points[mm - 1].z * dnorm) + zmin;
97 */
98 hz = h + zmin;
99 zz = points[mm - 1].z + zmin;
100 err = hz - zz;
101 xmm = points[mm - 1].x * dnorm + params->x_orig + west;
102 ymm = points[mm - 1].y * dnorm + params->y_orig + south;
103 if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig) &&
104 (ymm >= south + params->y_orig) && (ymm <= north + params->y_orig))
105 inside = 1;
106 else
107 inside = 0;
108
109 if (params->create_devi) {
110
111 if (inside) { /* if the point is inside the region */
114
116 cat = count;
117 Vect_cat_set(Cats2, 1, cat);
119
121 sprintf(buf, "insert into %s values ( %d ", ff->table, cat);
123
124 sprintf(buf, ", %f", err);
126 db_append_string(&sql2, ")");
127 G_debug(3, "IL_check_at_points_2d: %s", db_get_string(&sql2));
128
132 G_fatal_error("Cannot insert new row: %s",
134 }
135 count++;
136 }
137 }
138 (*ertot) += err * err;
139 }
140
141 /* cv stuff */
142 if (params->cv) {
143 h = b[0];
144 for (m = 1; m <= n_points - 1; m++) {
145 xx = points[m - 1].x - skip_point.x;
146 yy = points[m - 1].y - skip_point.y;
147 r2 = yy * yy + xx * xx;
148 if (r2 != 0.) {
149 /* rfsta2 = fstar2 * r2; */
150 r = r2;
151 h = h + b[m] * params->interp(r, params->fi);
152 }
153 }
154 hz = h + zmin;
155 zz = skip_point.z + zmin;
156 skip_err = hz - zz;
157 xmm = skip_point.x * dnorm + params->x_orig + west;
158 ymm = skip_point.y * dnorm + params->y_orig + south;
159
160 if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig) &&
161 (ymm >= south + params->y_orig) && (ymm <= north + params->y_orig))
162 inside = 1;
163 else
164 inside = 0;
165
166 if (inside) { /* if the point is inside the region */
169
171 cat = count;
172 Vect_cat_set(Cats2, 1, cat);
174
176 sprintf(buf, "insert into %s values ( %d ", ff->table, cat);
178
179 sprintf(buf, ", %f", skip_err);
181 db_append_string(&sql2, ")");
182 G_debug(3, "IL_check_at_points_2d: %s", db_get_string(&sql2));
183
187 G_fatal_error("Cannot insert new row: %s",
189 }
190 count++;
191 }
192 } /* cv */
193
194 return 1;
195}
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
double b
double r
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
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
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)