GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
interp2d.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  * bug fixes(mask) and modif. for variable smoothing Mitasova Jan 1997
22  *
23  */
24 
25 
26 #include <stdio.h>
27 #include <math.h>
28 #include <unistd.h>
29 #include <grass/gis.h>
30 #include <grass/glocale.h>
31 #include <grass/bitmap.h>
32 
33 #include <grass/interpf.h>
34 
35 
36 #define CEULER .57721566
37 
38 
39 int IL_grid_calc_2d(struct interp_params *params, struct quaddata *data, /* given segment */
40  struct BM *bitmask, /* bitmask */
41  double zmin, double zmax, /* min and max input z-values */
42  double *zminac, double *zmaxac, /* min and max interp. z-values */
43  double *gmin, double *gmax, /* min and max inperp. slope val. */
44  double *c1min, double *c1max, double *c2min, double *c2max, /* min and max interp. curv. val. */
45  double *ertot, /* total interplating func. error */
46  double *b, /* solutions of linear equations */
47  int offset1, /* offset for temp file writing */
48  double dnorm)
49 
50 /*
51  * Calculates grid for the given segment represented by data (contains
52  * n_rows, n_cols, ew_res,ns_res, and all points inside + overlap) using
53  * solutions of system of lin. equations and interpolating functions
54  * interp() and interpder(). Also calls secpar() to compute slope, aspect
55  * and curvatures if required.
56  */
57 {
58 
59  /*
60  * C C INTERPOLATION BY FUNCTIONAL METHOD : TPS + complete regul.
61  * c
62  */
63  double x_or = data->x_orig;
64  double y_or = data->y_orig;
65  int n_rows = data->n_rows;
66  int n_cols = data->n_cols;
67  int n_points = data->n_points;
68  struct triple *points;
69  static double *w2 = NULL;
70  static double *w = NULL;
71  int cond1, cond2;
72  double r;
73  double stepix, stepiy, xx, xg, yg, xx2;
74  double rfsta2, /* cons, cons1, */ wm, dx, dy, dxx, dyy, dxy, h, bmgd1,
75  bmgd2;
76  double r2, gd1, gd2; /* for interpder() */
77  int n1, k, l, m;
78  int ngstc, nszc, ngstr, nszr;
79  double zz;
80  int bmask = 1;
81  static int first_time_z = 1;
82  int offset, offset2;
83  double fstar2 = params->fi * params->fi / 4.;
84  double tfsta2, tfstad;
85  double ns_res, ew_res;
86  double rsin = 0, rcos = 0, teta, scale = 0; /*anisotropy parameters - added by JH 2002 */
87  double xxr, yyr;
88 
89  if (params->theta) {
90  teta = params->theta / 57.295779; /* deg to rad */
91  rsin = sin(teta);
92  rcos = cos(teta);
93  }
94  if (params->scalex)
95  scale = params->scalex;
96 
97  ns_res = (((struct quaddata *)(data))->ymax -
98  ((struct quaddata *)(data))->y_orig) / data->n_rows;
99  ew_res = (((struct quaddata *)(data))->xmax -
100  ((struct quaddata *)(data))->x_orig) / data->n_cols;
101 
102  /* tfsta2 = fstar2 * 2.; modified after removing normalization of z */
103  tfsta2 = (fstar2 * 2.) / dnorm;
104  tfstad = tfsta2 / dnorm;
105  points = data->points;
106 
107  /*
108  * normalization
109  */
110  stepix = ew_res / dnorm;
111  stepiy = ns_res / dnorm;
112 
113  cond2 = ((params->adxx != NULL) || (params->adyy != NULL) ||
114  (params->adxy != NULL));
115  cond1 = ((params->adx != NULL) || (params->ady != NULL) || cond2);
116 
117  if (!w) {
118  if (!(w = (double *)G_malloc(sizeof(double) * (params->KMAX2 + 9)))) {
119  G_warning(_("Out of memory"));
120  return -1;
121  }
122  }
123  if (!w2) {
124  if (!(w2 = (double *)G_malloc(sizeof(double) * (params->KMAX2 + 9)))) {
125  G_warning(_("Out of memory"));
126  return -1;
127  }
128  }
129  n1 = n_points + 1;
130  /*
131  * C C INTERPOLATION * MOST INNER LOOPS ! C
132  */
133  ngstc = (int)(x_or / ew_res + 0.5) + 1;
134  nszc = ngstc + n_cols - 1;
135  ngstr = (int)(y_or / ns_res + 0.5) + 1;
136  nszr = ngstr + n_rows - 1;
137 
138 
139  for (k = ngstr; k <= nszr; k++) {
140  offset = offset1 * (k - 1); /* rows offset */
141  yg = (k - ngstr) * stepiy + stepiy / 2.; /* fixed by J.H. in July 01 */
142  for (m = 1; m <= n_points; m++) {
143  wm = yg - points[m - 1].y;
144  w[m] = wm;
145  w2[m] = wm * wm;
146  }
147  for (l = ngstc; l <= nszc; l++) {
148  if (bitmask != NULL)
149  /* if(params->maskmap != NULL) PK Apr 03 MASK support */
150  bmask = BM_get(bitmask, l - 1, k - 1); /*fixed by helena jan 97 */
151  /* if(bmask==0 || bmask==-1) fprintf(stderr, "bmask=%d, at (%d,%d)\n", bmask, l, k); */
152  xg = (l - ngstc) * stepix + stepix / 2.; /*fixed by J.H. in July 01 */
153  dx = 0.;
154  dy = 0.;
155  dxx = 0.;
156  dyy = 0.;
157  dxy = 0.;
158  zz = 0.;
159  if (bmask == 1) { /* compute everything for area which is
160  * not masked out */
161  h = b[0];
162  for (m = 1; m <= n_points; m++) {
163  xx = xg - points[m - 1].x;
164  if ((params->theta) && (params->scalex)) {
165  /* we run anisotropy */
166  xxr = xx * rcos + w[m] * rsin;
167  yyr = w[m] * rcos - xx * rsin;
168  xx2 = xxr * xxr;
169  w2[m] = yyr * yyr;
170  r2 = scale * xx2 + w2[m];
171  r = r2;
172  rfsta2 = scale * xx2 + w2[m];
173  }
174  else {
175  xx2 = xx * xx;
176  r2 = xx2 + w2[m];
177  r = r2;
178  rfsta2 = xx2 + w2[m];
179  }
180 
181  h = h + b[m] * params->interp(r, params->fi);
182  if (cond1) {
183  if (!params->interpder(r, params->fi, &gd1, &gd2))
184  return -1;
185  bmgd1 = b[m] * gd1;
186  dx = dx + bmgd1 * xx;
187  dy = dy + bmgd1 * w[m];
188  if (cond2) {
189  bmgd2 = b[m] * gd2;
190  dxx = dxx + bmgd2 * xx2 + bmgd1;
191  dyy = dyy + bmgd2 * w2[m] + bmgd1;
192  dxy = dxy + bmgd2 * xx * w[m];
193  }
194  }
195  }
196 
197  /* zz = (h * dnorm) + zmin; replaced by helena jan. 97 due to removing norma
198  lization of z and zm in segmen2d.c */
199  zz = h + zmin;
200  if (first_time_z) {
201  first_time_z = 0;
202  *zmaxac = *zminac = zz;
203  }
204  *zmaxac = amax1(zz, *zmaxac);
205  *zminac = amin1(zz, *zminac);
206  if ((zz > zmax + 0.1 * (zmax - zmin))
207  || (zz < zmin - 0.1 * (zmax - zmin))) {
208  static int once = 0;
209 
210  if (!once) {
211  once = 1;
212  G_warning(_("Overshoot - increase in tension suggested. "
213  "Overshoot occures at (%d,%d) cell. "
214  "Z-value %f, zmin %f, zmax %f."),
215  l, k, zz, zmin, zmax);
216  }
217  }
218 
219  params->az[l] = (FCELL) zz;
220 
221  if (cond1) {
222  params->adx[l] = (FCELL) (-dx * tfsta2);
223  params->ady[l] = (FCELL) (-dy * tfsta2);
224  if (cond2) {
225  params->adxx[l] = (FCELL) (-dxx * tfstad);
226  params->adyy[l] = (FCELL) (-dyy * tfstad);
227  params->adxy[l] = (FCELL) (-dxy * tfstad);
228  }
229  }
230 
231  }
232  else {
233  G_set_d_null_value(params->az + l, 1);
234  /* fprintf (stderr, "zz=%f, az[l]=%f, c=%d\n", zz, params->az[l], l); */
235 
236  if (cond1) {
237  G_set_d_null_value(params->adx + l, 1);
238  G_set_d_null_value(params->ady + l, 1);
239  if (cond2) {
240  G_set_d_null_value(params->adxx + l, 1);
241  G_set_d_null_value(params->adyy + l, 1);
242  G_set_d_null_value(params->adxy + l, 1);
243  }
244  }
245  }
246 
247  }
248  if (cond1 && (params->deriv != 1)) {
249  if (params->secpar(params, ngstc, nszc, k, bitmask,
250  gmin, gmax, c1min, c1max, c2min, c2max, cond1,
251  cond2) < 0)
252  return -1;
253  }
254 
255  offset2 = (offset + ngstc - 1) * sizeof(FCELL);
256  if (params->wr_temp(params, ngstc, nszc, offset2) < 0)
257  return -1;
258  }
259  return 1;
260 }
int BM_get(struct BM *map, int x, int y)
Gets &#39;val&#39; from the bitmap.
Definition: bitmap.c:220
int l
Definition: dataquad.c:292
float b
Definition: named_colr.c:8
double y_orig
Definition: dataquad.h:33
int(* secpar)()
Definition: interpf.h:67
double scalex
Definition: interpf.h:58
void G_set_d_null_value(DCELL *dcellVals, int numVals)
Definition: null_val.c:176
float r
Definition: named_colr.c:8
int(* wr_temp)()
Definition: interpf.h:70
double theta
Definition: interpf.h:57
double amax1(double, double)
Definition: minmax.c:54
DCELL * adxx
Definition: interpf.h:48
int(* interpder)()
Definition: interpf.h:69
double(* interp)()
Definition: interpf.h:68
double x_orig
Definition: dataquad.h:32
DCELL * adx
Definition: interpf.h:48
struct triple * points
Definition: dataquad.h:39
tuple data
DCELL * az
Definition: interpf.h:48
double x
Definition: dataquad.h:24
int
Definition: g3dcolor.c:48
double fi
Definition: interpf.h:49
double ymax
Definition: dataquad.c:293
double amin1(double, double)
Definition: minmax.c:67
DCELL * ady
Definition: interpf.h:48
return NULL
Definition: dbfopen.c:1394
G_warning("category support for [%s] in mapset [%s] %s", name, mapset, type)
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, double *b, int offset1, double dnorm)
Definition: interp2d.c:39
int n_rows
Definition: dataquad.h:36
DCELL * adxy
Definition: interpf.h:48
int n_cols
Definition: dataquad.h:37
double y
Definition: dataquad.h:25
tuple h
panel.defaultSize = wx.CheckBox(panel, id = wx.ID_ANY, label = _(&quot;Use default size&quot;)) panel...
int n_points
Definition: dataquad.h:38
DCELL * adyy
Definition: interpf.h:48