GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
interp2d.c
Go to the documentation of this file.
1/*!
2 * \file interp2d.c
3 *
4 * \author
5 * Lubos Mitas (original program and various modifications)
6 *
7 * \author
8 * H. Mitasova,
9 * I. Kosinovsky, D. Gerdes,
10 * D. McCauley
11 * (GRASS4.1 version of the program and GRASS4.2 modifications)
12 *
13 * \author
14 * L. Mitas,
15 * H. Mitasova,
16 * I. Kosinovsky,
17 * D.Gerdes,
18 * D. McCauley
19 * (1993, 1995)
20 *
21 * \author modified by McCauley in August 1995
22 * \author modified by Mitasova in August 1995, Nov. 1996
23 * \author
24 * bug fixes(mask) and modification for variable smoothing
25 * Mitasova (Jan 1997)
26 *
27 * \copyright
28 * (C) 1993-1999 by Lubos Mitas and the GRASS Development Team
29 *
30 * \copyright
31 * This program is free software under the
32 * GNU General Public License (>=v2).
33 * Read the file COPYING that comes with GRASS
34 * for details.
35 */
36
37#include <stdio.h>
38#include <math.h>
39#include <unistd.h>
40
41#include <grass/gis.h>
42#include <grass/raster.h>
43#include <grass/glocale.h>
44#include <grass/bitmap.h>
45
46#include <grass/interpf.h>
47
48#define CEULER .57721566
49
50/*!
51 * Calculates grid values for a given segment
52 *
53 * Calculates grid for the given segment represented by data (contains
54 * n_rows, n_cols, ew_res,ns_res, and all points inside + overlap) using
55 * solutions of system of linear equations and interpolating functions
56 * interp() and interpder(). Also calls secpar() to compute slope, aspect
57 * and curvatures if required.
58 *
59 * *ertot* can be also called *RMS deviation of the interpolated surface*
60 */
62 struct interp_params *params, struct quaddata *data, /*!< given segment */
63 struct BM *bitmask, /*!< bitmask */
64 double zmin, double zmax, /*!< min and max input z-values */
65 double *zminac, double *zmaxac, /*!< min and max interp. z-values */
66 double *gmin, double *gmax, /*!< min and max interp. slope val. */
67 double *c1min, double *c1max, /*!< min and max interp. curv. val. */
68 double *c2min, double *c2max, /*!< min and max interp. curv. val. */
69 double *ertot UNUSED, /*!< total interpolating func. error */
70 double *b, /*!< solutions of linear equations */
71 off_t offset1, /*!< offset for temp file writing */
72 double dnorm)
73{
74
75 /*
76 * C C INTERPOLATION BY FUNCTIONAL METHOD : TPS + complete regul.
77 * c
78 */
79 double x_or = data->x_orig;
80 double y_or = data->y_orig;
81 int n_rows = data->n_rows;
82 int n_cols = data->n_cols;
83 int n_points = data->n_points;
84 struct triple *points;
85 static double *w2 = NULL;
86 static double *w = NULL;
87 int cond1, cond2;
88 double r;
89 double stepix, stepiy, xx, xg, yg, xx2;
90 double /* rfsta2, cons, cons1, */ wm, dx, dy, dxx, dyy, dxy, h, bmgd1,
91 bmgd2;
92 double r2, gd1, gd2; /* for interpder() */
93 int /* n1, */ k, l, m;
94 int ngstc, nszc, ngstr, nszr;
95 double zz;
96 int bmask = 1;
97 static int first_time_z = 1;
98 off_t offset, offset2;
99 double fstar2 = params->fi * params->fi / 4.;
100 double tfsta2, tfstad;
101 double ns_res, ew_res;
102 double rsin = 0, rcos = 0, teta,
103 scale = 0; /*anisotropy parameters - added by JH 2002 */
104 double xxr, yyr;
105
106 if (params->theta) {
107 teta = params->theta / M_R2D; /* deg to rad */
108 rsin = sin(teta);
109 rcos = cos(teta);
110 }
111 if (params->scalex)
112 scale = params->scalex;
113
114 ns_res = (((struct quaddata *)(data))->ymax -
115 ((struct quaddata *)(data))->y_orig) /
116 data->n_rows;
117 ew_res = (((struct quaddata *)(data))->xmax -
118 ((struct quaddata *)(data))->x_orig) /
119 data->n_cols;
120
121 /* tfsta2 = fstar2 * 2.; modified after removing normalization of z */
122 tfsta2 = (fstar2 * 2.) / dnorm;
123 tfstad = tfsta2 / dnorm;
124 points = data->points;
125
126 /*
127 * normalization
128 */
129 stepix = ew_res / dnorm;
130 stepiy = ns_res / dnorm;
131
132 cond2 = ((params->adxx != NULL) || (params->adyy != NULL) ||
133 (params->adxy != NULL));
134 cond1 = ((params->adx != NULL) || (params->ady != NULL) || cond2);
135
136 if (!w) {
137 if (!(w = (double *)G_malloc(sizeof(double) * (params->KMAX2 + 9)))) {
138 G_warning(_("Out of memory"));
139 return -1;
140 }
141 }
142 if (!w2) {
143 if (!(w2 = (double *)G_malloc(sizeof(double) * (params->KMAX2 + 9)))) {
144 G_warning(_("Out of memory"));
145 return -1;
146 }
147 }
148 /* n1 = n_points + 1; */
149 /*
150 * C C INTERPOLATION * MOST INNER LOOPS ! C
151 */
152 ngstc = (int)(x_or / ew_res + 0.5) + 1;
153 nszc = ngstc + n_cols - 1;
154 ngstr = (int)(y_or / ns_res + 0.5) + 1;
155 nszr = ngstr + n_rows - 1;
156
157 for (k = ngstr; k <= nszr; k++) {
158 offset = offset1 * (k - 1); /* rows offset */
159 yg = (k - ngstr) * stepiy + stepiy / 2.; /* fixed by J.H. in July 01 */
160 for (m = 1; m <= n_points; m++) {
161 wm = yg - points[m - 1].y;
162 w[m] = wm;
163 w2[m] = wm * wm;
164 }
165 for (l = ngstc; l <= nszc; l++) {
166 if (bitmask != NULL)
167 /* if(params->maskmap != NULL) PK Apr 03 MASK support */
168 bmask =
169 BM_get(bitmask, l - 1, k - 1); /*fixed by helena jan 97 */
170 /* if(bmask==0 || bmask==-1) fprintf(stderr, "bmask=%d, at
171 * (%d,%d)\n", bmask, l, k); */
172 xg = (l - ngstc) * stepix +
173 stepix / 2.; /*fixed by J.H. in July 01 */
174 dx = 0.;
175 dy = 0.;
176 dxx = 0.;
177 dyy = 0.;
178 dxy = 0.;
179 zz = 0.;
180 if (bmask == 1) { /* compute everything for area which is
181 * not masked out */
182 h = b[0];
183 for (m = 1; m <= n_points; m++) {
184 xx = xg - points[m - 1].x;
185 if ((params->theta) && (params->scalex)) {
186 /* we run anisotropy */
187 xxr = xx * rcos + w[m] * rsin;
188 yyr = w[m] * rcos - xx * rsin;
189 xx2 = xxr * xxr;
190 w2[m] = yyr * yyr;
191 r2 = scale * xx2 + w2[m];
192 r = r2;
193 /* rfsta2 = scale * xx2 + w2[m]; */
194 }
195 else {
196 xx2 = xx * xx;
197 r2 = xx2 + w2[m];
198 r = r2;
199 /* rfsta2 = xx2 + w2[m]; */
200 }
201
202 h = h + b[m] * params->interp(r, params->fi);
203 if (cond1) {
204 if (!params->interpder(r, params->fi, &gd1, &gd2))
205 return -1;
206 bmgd1 = b[m] * gd1;
207 dx = dx + bmgd1 * xx;
208 dy = dy + bmgd1 * w[m];
209 if (cond2) {
210 bmgd2 = b[m] * gd2;
211 dxx = dxx + bmgd2 * xx2 + bmgd1;
212 dyy = dyy + bmgd2 * w2[m] + bmgd1;
213 dxy = dxy + bmgd2 * xx * w[m];
214 }
215 }
216 }
217
218 /* zz = (h * dnorm) + zmin; replaced by helena jan. 97 due
219 to removing norma lization of z and zm in segmen2d.c */
220 zz = h + zmin;
221 if (first_time_z) {
222 first_time_z = 0;
223 *zmaxac = *zminac = zz;
224 }
225 *zmaxac = amax1(zz, *zmaxac);
226 *zminac = amin1(zz, *zminac);
227 if ((zz > zmax + 0.1 * (zmax - zmin)) ||
228 (zz < zmin - 0.1 * (zmax - zmin))) {
229 static int once = 0;
230
231 if (!once) {
232 once = 1;
233 G_warning(
234 _("Overshoot - increase in tension suggested. "
235 "Overshoot occurs at (%d,%d) cell. "
236 "Z-value %f, zmin %f, zmax %f."),
237 l, k, zz, zmin, zmax);
238 }
239 }
240
241 params->az[l] = (FCELL)zz;
242
243 if (cond1) {
244 params->adx[l] = (FCELL)(-dx * tfsta2);
245 params->ady[l] = (FCELL)(-dy * tfsta2);
246 if (cond2) {
247 params->adxx[l] = (FCELL)(-dxx * tfstad);
248 params->adyy[l] = (FCELL)(-dyy * tfstad);
249 params->adxy[l] = (FCELL)(-dxy * tfstad);
250 }
251 }
252 }
253 else {
254 Rast_set_d_null_value(params->az + l, 1);
255 /* fprintf (stderr, "zz=%f, az[l]=%f, c=%d\n", zz,
256 * params->az[l], l); */
257
258 if (cond1) {
259 Rast_set_d_null_value(params->adx + l, 1);
260 Rast_set_d_null_value(params->ady + l, 1);
261 if (cond2) {
262 Rast_set_d_null_value(params->adxx + l, 1);
263 Rast_set_d_null_value(params->adyy + l, 1);
264 Rast_set_d_null_value(params->adxy + l, 1);
265 }
266 }
267 }
268 }
269 if (cond1 && (params->deriv != 1)) {
270 if (params->secpar(params, ngstc, nszc, k, bitmask, gmin, gmax,
271 c1min, c1max, c2min, c2max, cond1, cond2) < 0)
272 return -1;
273 }
274
275 offset2 = (offset + ngstc - 1) * sizeof(FCELL);
276 if (params->wr_temp(params, ngstc, nszc, offset2) < 0)
277 return -1;
278 }
279 return 1;
280}
int BM_get(struct BM *map, int x, int y)
Gets 'val' from the bitmap.
Definition bitmap.c:217
#define NULL
Definition ccmath.h:32
double b
double l
double r
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
int IL_grid_calc_2d(struct interp_params *params, struct quaddata *data, struct BM *bitmask, double zmin, double zmax, double *zminac, double *zmaxac, double *gmin, double *gmax, double *c1min, double *c1max, double *c2min, double *c2max, double *ertot UNUSED, double *b, off_t offset1, double dnorm)
Definition interp2d.c:61
double amin1(double, double)
Definition minmax.c:65
double amax1(double, double)
Definition minmax.c:52
DCELL * az
Definition interpf.h:85
interp_fn * interp
Definition interpf.h:124
secpar_fn * secpar
Definition interpf.h:122
double fi
Definition interpf.h:88
double theta
Definition interpf.h:105
DCELL * adxy
Definition interpf.h:85
DCELL * adyy
Definition interpf.h:85
DCELL * adx
Definition interpf.h:85
DCELL * ady
Definition interpf.h:85
double scalex
Definition interpf.h:107
wr_temp_fn * wr_temp
Definition interpf.h:128
interpder_fn * interpder
Definition interpf.h:126
DCELL * adxx
Definition interpf.h:85
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
int n_cols
Definition dataquad.h:51
int n_rows
Definition dataquad.h:50
double x
Definition dataquad.h:39
double y
Definition dataquad.h:40