GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r73122
segmen2d.c
Go to the documentation of this file.
1 /*!
2  * \file segmen2d.c
3  *
4  * \author H. Mitasova, I. Kosinovsky, D. Gerdes
5  *
6  * \copyright
7  * (C) 1993 by Helena Mitasova and the GRASS Development Team
8  *
9  * \copyright
10  * This program is free software under the
11  * GNU General Public License (>=v2).
12  * Read the file COPYING that comes with GRASS
13  * for details.
14  *
15  */
16
17
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <math.h>
21 #include <grass/gis.h>
22 #include <grass/glocale.h>
23 #include <grass/interpf.h>
24 #include <grass/gmath.h>
25
26
27 /*!
28  * Interpolate recursively a tree of segments
29  *
30  * Recursively processes each segment in a tree by:
31  * - finding points from neighbouring segments so that the total number of
32  * points is between KMIN and KMAX2 by calling tree function MT_get_region().
33  * - creating and solving the system of linear equations using these points
34  * and interp() by calling matrix_create() and G_ludcmp().
35  * - checking the interpolating function values at points by calling
36  * check_points().
37  * - computing grid for this segment using points and interp() by calling
38  * grid_calc().
39  *
40  * \todo
41  * Isn't this in fact the updated version of the function (IL_interp_segments_new_2d)?
42  * The function IL_interp_segments_new_2d has the following, better behavior:
43  * The difference between this function and IL_interp_segments_2d() is making
44  * sure that additional points are taken from all directions, i.e. it finds
45  * equal number of points from neighboring segments in each of 8 neighborhoods.
46  */
48  struct tree_info *info, /*!< info for the quad tree */
49  struct multtree *tree, /*!< current leaf of the quad tree */
50  struct BM *bitmask, /*!< bitmask */
51  double zmin, double zmax, /*!< min and max input z-values */
52  double *zminac, double *zmaxac, /*!< min and max interp. z-values */
53  double *gmin, double *gmax, /*!< min and max inperp. slope val. */
54  double *c1min, double *c1max, /*!< min and max interp. curv. val. */
55  double *c2min, double *c2max, /*!< min and max interp. curv. val. */
56  double *ertot, /*!< total interplating func. error */
57  int totsegm, /*!< total number of segments */
58  off_t offset1, /*!< offset for temp file writing */
59  double dnorm)
60 {
61  double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1, temp2;
62  int i, npt, nptprev, MAXENC;
63  struct quaddata *data;
64  static int cursegm = 0;
65  static double *b = NULL;
66  static int *indx = NULL;
67  static double **matrix = NULL;
68  double ew_res, ns_res;
69  static int first_time = 1;
70  static double smseg;
71  int MINPTS;
72  double pr;
73  struct triple *point;
74  struct triple skip_point;
75  int m_skip, skip_index, j, k, segtest;
76  double xx, yy, zz;
77
78  /* find the size of the smallest segment once */
79  if (first_time) {
80  smseg = smallest_segment(info->root, 4);
81  first_time = 0;
82  }
83  ns_res = (((struct quaddata *)(info->root->data))->ymax -
84  ((struct quaddata *)(info->root->data))->y_orig) /
85  params->nsizr;
86  ew_res =
87  (((struct quaddata *)(info->root->data))->xmax -
88  ((struct quaddata *)(info->root->data))->x_orig) / params->nsizc;
89
90  if (tree == NULL)
91  return -1;
92  if (tree->data == NULL)
93  return -1;
94  if (((struct quaddata *)(tree->data))->points == NULL) {
95  for (i = 0; i < 4; i++) {
96  IL_interp_segments_2d(params, info, tree->leafs[i],
97  bitmask, zmin, zmax, zminac, zmaxac, gmin,
98  gmax, c1min, c1max, c2min, c2max, ertot,
99  totsegm, offset1, dnorm);
100  }
101  return 1;
102  }
103  else {
104  distx = (((struct quaddata *)(tree->data))->n_cols * ew_res) * 0.1;
105  disty = (((struct quaddata *)(tree->data))->n_rows * ns_res) * 0.1;
106  distxp = 0;
107  distyp = 0;
108  xmn = ((struct quaddata *)(tree->data))->x_orig;
109  xmx = ((struct quaddata *)(tree->data))->xmax;
110  ymn = ((struct quaddata *)(tree->data))->y_orig;
111  ymx = ((struct quaddata *)(tree->data))->ymax;
112  i = 0;
113  MAXENC = 0;
114  /* data is a window with zero points; some fields don't make sense in this case
115  so they are zero (like resolution,dimensions */
116  /* CHANGE */
117  /* Calcutaing kmin for surrent segment (depends on the size) */
118
119 /*****if (smseg <= 0.00001) MINPTS=params->kmin; else {} ***/
120  pr = pow(2., (xmx - xmn) / smseg - 1.);
121  MINPTS =
122  params->kmin * (pr / (1 + params->kmin * pr / params->KMAX2));
123  /* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf, smseg=%lf, DX=%lf \n", MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */
124
125  data =
126  (struct quaddata *)quad_data_new(xmn - distx, ymn - disty,
127  xmx + distx, ymx + disty, 0, 0,
128  0, params->KMAX2);
129  npt = MT_region_data(info, info->root, data, params->KMAX2, 4);
130
131  while ((npt < MINPTS) || (npt > params->KMAX2)) {
132  if (i >= 70) {
133  G_warning(_("Taking too long to find points for interpolation - "
134  "please change the region to area where your points are. "
135  "Continuing calculations..."));
136  break;
137  }
138  i++;
139  if (npt > params->KMAX2)
140  /* decrease window */
141  {
142  MAXENC = 1;
143  nptprev = npt;
144  temp1 = distxp;
145  distxp = distx;
146  distx = distxp - fabs(distx - temp1) * 0.5;
147  temp2 = distyp;
148  distyp = disty;
149  disty = distyp - fabs(disty - temp2) * 0.5;
150  /* decrease by 50% of a previous change in window */
151  }
152  else {
153  nptprev = npt;
154  temp1 = distyp;
155  distyp = disty;
156  temp2 = distxp;
157  distxp = distx;
158  if (MAXENC) {
159  disty = fabs(disty - temp1) * 0.5 + distyp;
160  distx = fabs(distx - temp2) * 0.5 + distxp;
161  }
162  else {
163  distx += distx;
164  disty += disty;
165  }
166  /* decrease by 50% of extra distance */
167  }
168  data->x_orig = xmn - distx; /* update window */
169  data->y_orig = ymn - disty;
170  data->xmax = xmx + distx;
171  data->ymax = ymx + disty;
172  data->n_points = 0;
173  npt = MT_region_data(info, info->root, data, params->KMAX2, 4);
174  }
175
176  if (totsegm != 0) {
177  G_percent(cursegm, totsegm, 1);
178  }
179  data->n_rows = ((struct quaddata *)(tree->data))->n_rows;
180  data->n_cols = ((struct quaddata *)(tree->data))->n_cols;
181
182  /* for printing out overlapping segments */
183  ((struct quaddata *)(tree->data))->x_orig = xmn - distx;
184  ((struct quaddata *)(tree->data))->y_orig = ymn - disty;
185  ((struct quaddata *)(tree->data))->xmax = xmx + distx;
186  ((struct quaddata *)(tree->data))->ymax = ymx + disty;
187
188  data->x_orig = xmn;
189  data->y_orig = ymn;
190  data->xmax = xmx;
191  data->ymax = ymx;
192
193  if (!matrix) {
194  if (!
195  (matrix =
196  G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
197  G_warning(_("Out of memory"));
198  return -1;
199  }
200  }
201  if (!indx) {
202  if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) {
203  G_warning(_("Out of memory"));
204  return -1;
205  }
206  }
207  if (!b) {
208  if (!(b = G_alloc_vector(params->KMAX2 + 3))) {
209  G_warning(_("Out of memory"));
210  return -1;
211  }
212  }
213  /* allocate memory for CV points only if cv is performed */
214  if (params->cv) {
215  if (!
216  (point =
217  (struct triple *)G_malloc(sizeof(struct triple) *
218  data->n_points))) {
219  G_warning(_("Out of memory"));
220  return -1;
221  }
222  }
223
224  /*normalize the data so that the side of average segment is about 1m */
225  /* put data_points into point only if CV is performed */
226
227  for (i = 0; i < data->n_points; i++) {
228  data->points[i].x = (data->points[i].x - data->x_orig) / dnorm;
229  data->points[i].y = (data->points[i].y - data->y_orig) / dnorm;
230  if (params->cv) {
231  point[i].x = data->points[i].x; /*cv stuff */
232  point[i].y = data->points[i].y; /*cv stuff */
233  point[i].z = data->points[i].z; /*cv stuff */
234  }
235
236  /* commented out by Helena january 1997 as this is not necessary
237  although it may be useful to put normalization of z back?
238  data->points[i].z = data->points[i].z / dnorm;
239  this made smoothing self-adjusting based on dnorm
240  if (params->rsm < 0.) data->points[i].sm = data->points[i].sm / dnorm;
241  */
242  }
243
244  /* cv stuff */
245  if (params->cv)
246  m_skip = data->n_points;
247  else
248  m_skip = 1;
249
250  /* remove after cleanup - this is just for testing */
251  skip_point.x = 0.;
252  skip_point.y = 0.;
253  skip_point.z = 0.;
254
255
256  /*** TODO: parallelize this loop instead of the LU solver! ***/
257  for (skip_index = 0; skip_index < m_skip; skip_index++) {
258  if (params->cv) {
259  segtest = 0;
260  j = 0;
261  xx = point[skip_index].x * dnorm + data->x_orig +
262  params->x_orig;
263  yy = point[skip_index].y * dnorm + data->y_orig +
264  params->y_orig;
265  zz = point[skip_index].z;
266  if (xx >= data->x_orig + params->x_orig &&
267  xx <= data->xmax + params->x_orig &&
268  yy >= data->y_orig + params->y_orig &&
269  yy <= data->ymax + params->y_orig) {
270  segtest = 1;
271  skip_point.x = point[skip_index].x;
272  skip_point.y = point[skip_index].y;
273  skip_point.z = point[skip_index].z;
274  for (k = 0; k < m_skip; k++) {
275  if (k != skip_index && params->cv) {
276  data->points[j].x = point[k].x;
277  data->points[j].y = point[k].y;
278  data->points[j].z = point[k].z;
279  j++;
280  }
281  }
282  } /* segment area test */
283  }
284  if (!params->cv) {
285  if (params->
286  matrix_create(params, data->points, data->n_points,
287  matrix, indx) < 0)
288  return -1;
289  }
290  else if (segtest == 1) {
291  if (params->
292  matrix_create(params, data->points, data->n_points - 1,
293  matrix, indx) < 0)
294  return -1;
295  }
296  if (!params->cv) {
297  for (i = 0; i < data->n_points; i++)
298  b[i + 1] = data->points[i].z;
299  b[0] = 0.;
300  G_lubksb(matrix, data->n_points + 1, indx, b);
301  /* put here condition to skip error if not needed */
302  params->check_points(params, data, b, ertot, zmin, dnorm,
303  skip_point);
304  }
305  else if (segtest == 1) {
306  for (i = 0; i < data->n_points - 1; i++)
307  b[i + 1] = data->points[i].z;
308  b[0] = 0.;
309  G_lubksb(matrix, data->n_points, indx, b);
310  params->check_points(params, data, b, ertot, zmin, dnorm,
311  skip_point);
312  }
313  } /*end of cv loop */
314
315  if (!params->cv)
316  if ((params->Tmp_fd_z != NULL) || (params->Tmp_fd_dx != NULL) ||
317  (params->Tmp_fd_dy != NULL) || (params->Tmp_fd_xx != NULL) ||
318  (params->Tmp_fd_yy != NULL) || (params->Tmp_fd_xy != NULL)) {
319
320  if (params->grid_calc(params, data, bitmask,
321  zmin, zmax, zminac, zmaxac, gmin, gmax,
322  c1min, c1max, c2min, c2max, ertot, b,
323  offset1, dnorm) < 0)
324  return -1;
325  }
326
327  /* show after to catch 100% */
328  cursegm++;
329  if (totsegm < cursegm)
330  G_debug(1, "%d %d", totsegm, cursegm);
331
332  if (totsegm != 0) {
333  G_percent(cursegm, totsegm, 1);
334  }
335  /*
336  G_free_matrix(matrix);
337  G_free_ivector(indx);
338  G_free_vector(b);
339  */
340  G_free(data->points);
341  G_free(data);
342  }
343  return 1;
344 }
345
346 double smallest_segment(struct multtree *tree, int n_leafs)
347 {
348  static int first_time = 1;
349  int ii;
350  static double minside;
351  double side;
352
353  if (tree == NULL)
354  return 0;
355  if (tree->data == NULL)
356  return 0;
357  if (tree->leafs != NULL) {
358  for (ii = 0; ii < n_leafs; ii++) {
359  side = smallest_segment(tree->leafs[ii], n_leafs);
360  if (first_time) {
361  minside = side;
362  first_time = 0;
363  }
364  if (side < minside)
365  minside = side;
366  }
367  }
368  else {
369  side = ((struct quaddata *)(tree->data))->xmax -
370  ((struct quaddata *)(tree->data))->x_orig;
371  return side;
372  }
373
374  return minside;
375 }
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:149
double y_orig
Definition: dataquad.h:51
double smallest_segment(struct multtree *, int)
Definition: segmen2d.c:346
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition: dalloc.c:60
FILE * Tmp_fd_xy
Definition: interpf.h:92
FILE * Tmp_fd_yy
Definition: interpf.h:92
double z
Definition: dataquad.h:44
Definition: bitmap.h:17
double y_orig
Definition: interpf.h:87
#define NULL
Definition: ccmath.h:32
FILE * Tmp_fd_xx
Definition: interpf.h:92
double x_orig
Definition: dataquad.h:50
struct triple * points
Definition: dataquad.h:57
struct quaddata * data
Definition: qtree.h:58
int IL_interp_segments_2d(struct interp_params *, struct tree_info *, struct multtree *, struct BM *, double, double, double *, double *, double *, double *, double *, double *, double *, double *, double *, int, off_t, double)
Definition: segmen2d.c:47
struct quaddata * quad_data_new(double x_or, double y_or, double xmax, double ymax, int rows, int cols, int n_points, int kmax)
Definition: dataquad.c:62
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition: lu.c:104
grid_calc_fn * grid_calc
Definition: interpf.h:97
FILE * Tmp_fd_dx
Definition: interpf.h:92
double b
Definition: r_raster.c:39
check_points_fn * check_points
Definition: interpf.h:99
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
double x
Definition: dataquad.h:42
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition: percent.c:62
double x_orig
Definition: interpf.h:87
double ymax
Definition: dataquad.h:53
FILE * Tmp_fd_dy
Definition: interpf.h:92
FILE * Tmp_fd_z
Definition: interpf.h:92
int MT_region_data(struct tree_info *info, struct multtree *tree, struct quaddata *data, int MAX, int n_leafs)
Definition: qtree.c:194
int * G_alloc_ivector(size_t n)
Vector matrix memory allocation.
Definition: ialloc.c:41
#define _(str)
Definition: glocale.h:13
Definition: qtree.h:56
double xmax
Definition: dataquad.h:52
int n_rows
Definition: dataquad.h:54
struct multtree * root
Definition: qtree.h:53
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41
int n_cols
Definition: dataquad.h:55
double y
Definition: dataquad.h:43
struct multtree ** leafs
Definition: qtree.h:59
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
int n_points
Definition: dataquad.h:56