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