GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-bea8435a9e
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 }
#define NULL
Definition: ccmath.h:32
int BM_get(struct BM *, int, int)
Gets 'val' from the bitmap.
Definition: bitmap.c:217
void G_warning(const char *,...) __attribute__((format(printf
#define G_malloc(n)
Definition: defs/gis.h:94
void Rast_set_d_null_value(DCELL *, int)
To set a number of DCELL raster values to NULL.
Definition: null_val.c:153
float FCELL
Definition: gis.h:630
#define M_R2D
Definition: gis.h:167
#define UNUSED
A macro for an attribute, if attached to a variable, indicating that the variable is not used.
Definition: gis.h:47
#define _(str)
Definition: glocale.h:10
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
double b
Definition: r_raster.c:39
double l
Definition: r_raster.c:39
double r
Definition: r_raster.c:39
Definition: bitmap.h:17
DCELL * az
Definition: interpf.h:89
interp_fn * interp
Definition: interpf.h:131
secpar_fn * secpar
Definition: interpf.h:129
double fi
Definition: interpf.h:92
double theta
Definition: interpf.h:111
DCELL * adxy
Definition: interpf.h:89
DCELL * adyy
Definition: interpf.h:89
DCELL * adx
Definition: interpf.h:89
DCELL * ady
Definition: interpf.h:89
double scalex
Definition: interpf.h:114
wr_temp_fn * wr_temp
Definition: interpf.h:135
interpder_fn * interpder
Definition: interpf.h:133
DCELL * adxx
Definition: interpf.h:89
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
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