GRASS Programmer's Manual  6.4.3(2013)-r
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Macros Pages
point2d.c
Go to the documentation of this file.
1 
2 /*-
3  *
4  * Original program and various modifications:
5  * Lubos Mitas
6  *
7  * GRASS4.1 version of the program and GRASS4.2 modifications:
8  * H. Mitasova
9  * I. Kosinovsky, D. Gerdes
10  * D. McCauley
11  *
12  * Copyright 1993, 1995:
13  * L. Mitas ,
14  * H. Mitasova ,
15  * I. Kosinovsky, ,
16  * D.Gerdes
17  * D. McCauley
18  *
19  * modified by McCauley in August 1995
20  * modified by Mitasova in August 1995, Nov. 1996
21  *
22  */
23 
24 
25 #include <stdio.h>
26 #include <math.h>
27 #include <unistd.h>
28 #include <grass/gis.h>
29 #include <grass/site.h>
30 #include <grass/Vect.h>
31 #include <grass/dbmi.h>
32 
33 #define POINT2D_C
34 #include <grass/interpf.h>
35 
36 int IL_check_at_points_2d(struct interp_params *params, struct quaddata *data, /* current region */
37  double *b, /* solution of linear equations */
38  double *ertot, /* total error */
39  double zmin, /* min z-value */
40  double dnorm, struct triple skip_point)
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 {
48  int n_points = data->n_points; /* number of points */
49  struct triple *points = data->points; /* points for interpolation */
50  double east = data->xmax;
51  double west = data->x_orig;
52  double north = data->ymax;
53  double south = data->y_orig;
54  double rfsta2, errmax, h, xx, yy, r2, hz, zz, err, xmm, ymm, r;
55  double skip_err;
56  int n1, mm, m, cat;
57  double fstar2;
58  int inside;
59 
60  /* Site *site; */
61  char buf[1024];
62 
63 
64  /* if ((site = G_site_new_struct (-1, 2, 0, 1)) == NULL)
65  G_fatal_error ("Memory error for site struct"); */
66 
67  fstar2 = params->fi * params->fi / 4.;
68  errmax = .0;
69  n1 = n_points + 1;
70  for (mm = 1; mm <= n_points; mm++) {
71  h = b[0];
72  for (m = 1; m <= n_points; m++) {
73  xx = points[mm - 1].x - points[m - 1].x;
74  yy = points[mm - 1].y - points[m - 1].y;
75  r2 = yy * yy + xx * xx;
76  if (r2 != 0.) {
77  rfsta2 = fstar2 * r2;
78  r = r2;
79  h = h + b[m] * params->interp(r, params->fi);
80  }
81  }
82  /* modified by helena january 1997 - normalization of z was
83  removed from segm2d.c and interp2d.c
84  hz = (h * dnorm) + zmin;
85  zz = (points[mm - 1].z * dnorm) + zmin;
86  */
87  hz = h + zmin;
88  zz = points[mm - 1].z + zmin;
89  err = hz - zz;
90  xmm = points[mm - 1].x * dnorm + params->x_orig + west;
91  ymm = points[mm - 1].y * dnorm + params->y_orig + south;
92  if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig)
93  && (ymm >= south + params->y_orig) &&
94  (ymm <= north + params->y_orig))
95  inside = 1;
96  else
97  inside = 0;
98 
99  if (params->fddevi != NULL) {
100 
101  if (inside) { /* if the point is inside the region */
104 
105  Vect_append_point(Pnts, xmm, ymm, zz);
106  cat = count;
107  Vect_cat_set(Cats2, 1, cat);
108  Vect_write_line(&Map2, GV_POINT, Pnts, Cats2);
109 
111  sprintf(buf, "insert into %s values ( %d ", ff->table, cat);
112  db_append_string(&sql2, buf);
113 
114  sprintf(buf, ", %f", err);
115  db_append_string(&sql2, buf);
116  db_append_string(&sql2, ")");
117  G_debug(3, db_get_string(&sql2));
118 
119  if (db_execute_immediate(driver2, &sql2) != DB_OK) {
122  G_fatal_error("Cannot insert new row: %s",
123  db_get_string(&sql2));
124  }
125  count++;
126 
127  }
128  }
129  (*ertot) += err * err;
130  }
131 
132  /* cv stuff */
133  if (params->cv) {
134  h = b[0];
135  for (m = 1; m <= n_points - 1; m++) {
136  xx = points[m - 1].x - skip_point.x;
137  yy = points[m - 1].y - skip_point.y;
138  r2 = yy * yy + xx * xx;
139  if (r2 != 0.) {
140  rfsta2 = fstar2 * r2;
141  r = r2;
142  h = h + b[m] * params->interp(r, params->fi);
143  }
144  }
145  hz = h + zmin;
146  zz = skip_point.z + zmin;
147  skip_err = hz - zz;
148  xmm = skip_point.x * dnorm + params->x_orig + west;
149  ymm = skip_point.y * dnorm + params->y_orig + south;
150 
151  if ((xmm >= west + params->x_orig) && (xmm <= east + params->x_orig)
152  && (ymm >= south + params->y_orig) &&
153  (ymm <= north + params->y_orig))
154  inside = 1;
155  else
156  inside = 0;
157 
158  if (inside) { /* if the point is inside the region */
161 
162  Vect_append_point(Pnts, xmm, ymm, zz);
163  cat = count;
164  Vect_cat_set(Cats2, 1, cat);
165  Vect_write_line(&Map2, GV_POINT, Pnts, Cats2);
166 
168  sprintf(buf, "insert into %s values ( %d ", ff->table, cat);
169  db_append_string(&sql2, buf);
170 
171  sprintf(buf, ", %f", skip_err);
172  db_append_string(&sql2, buf);
173  db_append_string(&sql2, ")");
174  G_debug(3, db_get_string(&sql2));
175 
176  if (db_execute_immediate(driver2, &sql2) != DB_OK) {
179  G_fatal_error("Cannot insert new row: %s",
180  db_get_string(&sql2));
181  }
182  count++;
183  }
184  } /* cv */
185 
186 
187  return 1;
188 }