GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-112dd97adf
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 #include <stdio.h>
18 #include <stdlib.h>
19 #include <math.h>
20 #include <grass/gis.h>
21 #include <grass/glocale.h>
22 #include <grass/interpf.h>
23 #include <grass/gmath.h>
24 
25 /*!
26  * Interpolate recursively a tree of segments
27  *
28  * Recursively processes each segment in a tree by:
29  * - finding points from neighbouring segments so that the total number of
30  * points is between KMIN and KMAX2 by calling tree function MT_get_region().
31  * - creating and solving the system of linear equations using these points
32  * and interp() by calling matrix_create() and G_ludcmp().
33  * - checking the interpolating function values at points by calling
34  * check_points().
35  * - computing grid for this segment using points and interp() by calling
36  * grid_calc().
37  *
38  * \todo
39  * Isn't this in fact the updated version of the function
40  * (IL_interp_segments_new_2d)? The function IL_interp_segments_new_2d has the
41  * following, better behavior: The difference between this function and
42  * IL_interp_segments_2d() is making sure that additional points are taken from
43  * all directions, i.e. it finds equal number of points from neighboring
44  * segments in each of 8 neighborhoods.
45  */
47  struct interp_params *params,
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, 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 = (((struct quaddata *)(info->root->data))->xmax -
87  ((struct quaddata *)(info->root->data))->x_orig) /
88  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], bitmask, zmin,
97  zmax, zminac, zmaxac, gmin, gmax, c1min,
98  c1max, c2min, c2max, ertot, totsegm, offset1,
99  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
115  this case 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 = params->kmin * (pr / (1 + params->kmin * pr / params->KMAX2));
122  /* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf, smseg=%lf,
123  * DX=%lf \n", MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */
124 
125  data = (struct quaddata *)quad_data_new(xmn - distx, ymn - disty,
126  xmx + distx, ymx + disty, 0, 0,
127  0, params->KMAX2);
128  npt = MT_region_data(info, info->root, data, params->KMAX2, 4);
129 
130  while ((npt < MINPTS) || (npt > params->KMAX2)) {
131  if (i >= 70) {
132  G_warning(
133  _("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  temp1 = distxp;
144  distxp = distx;
145  distx = distxp - fabs(distx - temp1) * 0.5;
146  temp2 = distyp;
147  distyp = disty;
148  disty = distyp - fabs(disty - temp2) * 0.5;
149  /* decrease by 50% of a previous change in window */
150  }
151  else {
152  temp1 = distyp;
153  distyp = disty;
154  temp2 = distxp;
155  distxp = distx;
156  if (MAXENC) {
157  disty = fabs(disty - temp1) * 0.5 + distyp;
158  distx = fabs(distx - temp2) * 0.5 + distxp;
159  }
160  else {
161  distx += distx;
162  disty += disty;
163  }
164  /* decrease by 50% of extra distance */
165  }
166  data->x_orig = xmn - distx; /* update window */
167  data->y_orig = ymn - disty;
168  data->xmax = xmx + distx;
169  data->ymax = ymx + disty;
170  data->n_points = 0;
171  npt = MT_region_data(info, info->root, data, params->KMAX2, 4);
172  }
173 
174  if (totsegm != 0) {
175  G_percent(cursegm, totsegm, 1);
176  }
177  data->n_rows = ((struct quaddata *)(tree->data))->n_rows;
178  data->n_cols = ((struct quaddata *)(tree->data))->n_cols;
179 
180  /* for printing out overlapping segments */
181  ((struct quaddata *)(tree->data))->x_orig = xmn - distx;
182  ((struct quaddata *)(tree->data))->y_orig = ymn - disty;
183  ((struct quaddata *)(tree->data))->xmax = xmx + distx;
184  ((struct quaddata *)(tree->data))->ymax = ymx + disty;
185 
186  data->x_orig = xmn;
187  data->y_orig = ymn;
188  data->xmax = xmx;
189  data->ymax = ymx;
190 
191  if (!matrix) {
192  if (!(matrix =
193  G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
194  G_warning(_("Out of memory"));
195  return -1;
196  }
197  }
198  if (!indx) {
199  if (!(indx = G_alloc_ivector(params->KMAX2 + 1))) {
200  G_warning(_("Out of memory"));
201  return -1;
202  }
203  }
204  if (!b) {
205  if (!(b = G_alloc_vector(params->KMAX2 + 3))) {
206  G_warning(_("Out of memory"));
207  return -1;
208  }
209  }
210  /* allocate memory for CV points only if cv is performed */
211  if (params->cv) {
212  if (!(point = (struct triple *)G_malloc(sizeof(struct triple) *
213  data->n_points))) {
214  G_warning(_("Out of memory"));
215  return -1;
216  }
217  }
218 
219  /*normalize the data so that the side of average segment is about 1m */
220  /* put data_points into point only if CV is performed */
221 
222  for (i = 0; i < data->n_points; i++) {
223  data->points[i].x = (data->points[i].x - data->x_orig) / dnorm;
224  data->points[i].y = (data->points[i].y - data->y_orig) / dnorm;
225  if (params->cv) {
226  point[i].x = data->points[i].x; /*cv stuff */
227  point[i].y = data->points[i].y; /*cv stuff */
228  point[i].z = data->points[i].z; /*cv stuff */
229  }
230 
231  /* commented out by Helena january 1997 as this is not necessary
232  although it may be useful to put normalization of z back?
233  data->points[i].z = data->points[i].z / dnorm;
234  this made smoothing self-adjusting based on dnorm
235  if (params->rsm < 0.) data->points[i].sm = data->points[i].sm /
236  dnorm;
237  */
238  }
239 
240  /* cv stuff */
241  if (params->cv)
242  m_skip = data->n_points;
243  else
244  m_skip = 1;
245 
246  /* remove after cleanup - this is just for testing */
247  skip_point.x = 0.;
248  skip_point.y = 0.;
249  skip_point.z = 0.;
250 
251  /*** TODO: parallelize this loop instead of the LU solver! ***/
252  for (skip_index = 0; skip_index < m_skip; skip_index++) {
253  if (params->cv) {
254  segtest = 0;
255  j = 0;
256  xx =
257  point[skip_index].x * dnorm + data->x_orig + params->x_orig;
258  yy =
259  point[skip_index].y * dnorm + data->y_orig + params->y_orig;
260  /* zz = point[skip_index].z; */
261  if (xx >= data->x_orig + params->x_orig &&
262  xx <= data->xmax + params->x_orig &&
263  yy >= data->y_orig + params->y_orig &&
264  yy <= data->ymax + params->y_orig) {
265  segtest = 1;
266  skip_point.x = point[skip_index].x;
267  skip_point.y = point[skip_index].y;
268  skip_point.z = point[skip_index].z;
269  for (k = 0; k < m_skip; k++) {
270  if (k != skip_index && params->cv) {
271  data->points[j].x = point[k].x;
272  data->points[j].y = point[k].y;
273  data->points[j].z = point[k].z;
274  j++;
275  }
276  }
277  } /* segment area test */
278  }
279  if (!params->cv) {
280  if (params->matrix_create(params, data->points, data->n_points,
281  matrix, indx) < 0)
282  return -1;
283  }
284  else if (segtest == 1) {
285  if (params->matrix_create(params, data->points,
286  data->n_points - 1, matrix, indx) < 0)
287  return -1;
288  }
289  if (!params->cv) {
290  for (i = 0; i < data->n_points; i++)
291  b[i + 1] = data->points[i].z;
292  b[0] = 0.;
293  G_lubksb(matrix, data->n_points + 1, indx, b);
294  /* put here condition to skip error if not needed */
295  params->check_points(params, data, b, ertot, zmin, dnorm,
296  skip_point);
297  }
298  else if (segtest == 1) {
299  for (i = 0; i < data->n_points - 1; i++)
300  b[i + 1] = data->points[i].z;
301  b[0] = 0.;
302  G_lubksb(matrix, data->n_points, indx, b);
303  params->check_points(params, data, b, ertot, zmin, dnorm,
304  skip_point);
305  }
306  } /*end of cv loop */
307 
308  if (!params->cv)
309  if ((params->Tmp_fd_z != NULL) || (params->Tmp_fd_dx != NULL) ||
310  (params->Tmp_fd_dy != NULL) || (params->Tmp_fd_xx != NULL) ||
311  (params->Tmp_fd_yy != NULL) || (params->Tmp_fd_xy != NULL)) {
312 
313  if (params->grid_calc(params, data, bitmask, zmin, zmax, zminac,
314  zmaxac, gmin, gmax, c1min, c1max, c2min,
315  c2max, ertot, b, offset1, dnorm) < 0)
316  return -1;
317  }
318 
319  /* show after to catch 100% */
320  cursegm++;
321  if (totsegm < cursegm)
322  G_debug(1, "%d %d", totsegm, cursegm);
323 
324  if (totsegm != 0) {
325  G_percent(cursegm, totsegm, 1);
326  }
327  /*
328  G_free_matrix(matrix);
329  G_free_ivector(indx);
330  G_free_vector(b);
331  */
332  G_free(data->points);
333  G_free(data);
334  }
335  return 1;
336 }
337 
338 double smallest_segment(struct multtree *tree, int n_leafs)
339 {
340  static int first_time = 1;
341  int ii;
342  static double minside;
343  double side;
344 
345  if (tree == NULL)
346  return 0;
347  if (tree->data == NULL)
348  return 0;
349  if (tree->leafs != NULL) {
350  for (ii = 0; ii < n_leafs; ii++) {
351  side = smallest_segment(tree->leafs[ii], n_leafs);
352  if (first_time) {
353  minside = side;
354  first_time = 0;
355  }
356  if (side < minside)
357  minside = side;
358  }
359  }
360  else {
361  side = ((struct quaddata *)(tree->data))->xmax -
362  ((struct quaddata *)(tree->data))->x_orig;
363  return side;
364  }
365 
366  return minside;
367 }
#define NULL
Definition: ccmath.h:32
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:59
void G_percent(long, long, int)
Print percent complete messages.
Definition: percent.c:61
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
void G_warning(const char *,...) __attribute__((format(printf
#define G_malloc(n)
Definition: defs/gis.h:94
int G_debug(int, const char *,...) __attribute__((format(printf
int * G_alloc_ivector(size_t)
Vector matrix memory allocation.
Definition: ialloc.c:38
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition: lu.c:103
double * G_alloc_vector(size_t)
Vector matrix memory allocation.
Definition: dalloc.c:39
double ** G_alloc_matrix(int, int)
Matrix memory allocation.
Definition: dalloc.c:57
#define _(str)
Definition: glocale.h:10
int MT_region_data(struct tree_info *info, struct multtree *tree, struct quaddata *data, int MAX, int n_leafs)
Definition: qtree.c:186
double b
Definition: r_raster.c:39
double smallest_segment(struct multtree *tree, int n_leafs)
Definition: segmen2d.c:338
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:46
Definition: bitmap.h:17
check_points_fn * check_points
Definition: interpf.h:120
FILE * Tmp_fd_xx
Definition: interpf.h:111
FILE * Tmp_fd_xy
Definition: interpf.h:111
FILE * Tmp_fd_yy
Definition: interpf.h:111
grid_calc_fn * grid_calc
Definition: interpf.h:116
double x_orig
Definition: interpf.h:101
FILE * Tmp_fd_dx
Definition: interpf.h:111
double y_orig
Definition: interpf.h:101
FILE * Tmp_fd_z
Definition: interpf.h:111
FILE * Tmp_fd_dy
Definition: interpf.h:111
matrix_create_fn * matrix_create
Definition: interpf.h:118
Definition: qtree.h:52
struct multtree ** leafs
Definition: qtree.h:54
struct quaddata * data
Definition: qtree.h:53
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
struct multtree * root
Definition: qtree.h:49
double z
Definition: dataquad.h:41
double x
Definition: dataquad.h:39
double y
Definition: dataquad.h:40