GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-cb7cbab7e7
segmen2d_parallel.c
Go to the documentation of this file.
1 /*!
2  * \file segmen2d.c
3  *
4  * \author H. Mitasova, I. Kosinovsky, D. Gerdes (single core)
5  * \author Stanislav Zubal, Michal Lacko (OpenMP version)
6  * \author Anna Petrasova (OpenMP version GRASS integration)
7  *
8  * \copyright
9  * (C) 1993-2017 by Helena Mitasova and the GRASS Development Team
10  *
11  * \copyright
12  * This program is free software under the
13  * GNU General Public License (>=v2).
14  * Read the file COPYING that comes with GRASS
15  * for details.
16  *
17  */
18 
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <math.h>
22 #if defined(_OPENMP)
23 #include <omp.h>
24 #endif
25 #include <grass/gis.h>
26 #include <grass/glocale.h>
27 #include <grass/interpf.h>
28 #include <grass/gmath.h>
29 
30 static int cut_tree(struct multtree *, struct multtree **, int *);
31 
32 /*!
33  * See documentation for IL_interp_segments_2d.
34  * This is a parallel processing implementation.
35  */
37  struct interp_params *params,
38  struct tree_info *info, /*!< info for the quad tree */
39  struct multtree *tree, /*!< current leaf of the quad tree */
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, /*!< min and max interp. curv. val. */
45  double *c2min, double *c2max, /*!< min and max interp. curv. val. */
46  double *ertot, /*!< total interplating func. error */
47  int totsegm, /*!< total number of segments */
48  off_t offset1, /*!< offset for temp file writing */
49  double dnorm, int threads)
50 {
51  int some_thread_failed = 0;
52  int tid = 0;
53  int i = 0;
54  int j = 0;
55  int i_cnt;
56  int cursegm = 0;
57  double smseg;
58  double ***matrix = NULL;
59  int **indx = NULL;
60  double **b = NULL;
61  double **A = NULL;
62  struct quaddata **data_local;
63  struct multtree **all_leafs;
64 
65  all_leafs =
66  (struct multtree **)G_malloc(sizeof(struct multtree *) * totsegm);
67  data_local =
68  (struct quaddata **)G_malloc(sizeof(struct quaddata *) * threads);
69  matrix = (double ***)G_malloc(sizeof(double **) * threads);
70  indx = (int **)G_malloc(sizeof(int *) * threads);
71  b = (double **)G_malloc(sizeof(double *) * threads);
72  A = (double **)G_malloc(sizeof(double *) * threads);
73 
74  for (i_cnt = 0; i_cnt < threads; i_cnt++) {
75  if (!(matrix[i_cnt] =
76  G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
77  G_fatal_error(_("Out of memory"));
78  return -1;
79  }
80  }
81 
82  for (i_cnt = 0; i_cnt < threads; i_cnt++) {
83  if (!(indx[i_cnt] = G_alloc_ivector(params->KMAX2 + 1))) {
84  G_fatal_error(_("Out of memory"));
85  return -1;
86  }
87  }
88 
89  for (i_cnt = 0; i_cnt < threads; i_cnt++) {
90  if (!(b[i_cnt] = G_alloc_vector(params->KMAX2 + 3))) {
91  G_fatal_error(_("Out of memory"));
92  return -1;
93  }
94  }
95 
96  for (i_cnt = 0; i_cnt < threads; i_cnt++) {
97  if (!(A[i_cnt] = G_alloc_vector(
98  (params->KMAX2 + 2) * (params->KMAX2 + 2) + 1))) {
99  G_fatal_error(_("Out of memory"));
100  return -1;
101  }
102  }
103 
104  smseg = smallest_segment(tree, 4);
105  cut_tree(tree, all_leafs, &i);
106 
107  G_message(_("Starting parallel work"));
108 #pragma omp parallel firstprivate( \
109  tid, i, j, zmin, zmax, tree, totsegm, offset1, dnorm, smseg, ertot, \
110  params, info, all_leafs, bitmask, b, indx, matrix, data_local, A) \
111  shared(cursegm, threads, some_thread_failed, zminac, zmaxac, gmin, gmax, \
112  c1min, c1max, c2min, c2max) default(none)
113  {
114 #pragma omp for schedule(dynamic)
115  for (i_cnt = 0; i_cnt < totsegm; i_cnt++) {
116  /* Obtain thread id */
117 #if defined(_OPENMP)
118  tid = omp_get_thread_num();
119 #endif
120 
121  double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1,
122  temp2;
123  int npt, MAXENC;
124  double ew_res, ns_res;
125  int MINPTS;
126  double pr;
127  struct triple *point;
128  struct triple target_point;
129  int npoints, point_index, k;
130  double xx, yy /*, zz */;
131  double xmm, ymm, err, pointz;
132 
133  // struct quaddata *data_local;
134 
135  ns_res = (((struct quaddata *)(tree->data))->ymax -
136  ((struct quaddata *)(tree->data))->y_orig) /
137  params->nsizr;
138  ew_res = (((struct quaddata *)(tree->data))->xmax -
139  ((struct quaddata *)(tree->data))->x_orig) /
140  params->nsizc;
141 
142  if (all_leafs[i_cnt] == NULL) {
143  some_thread_failed = -1;
144  continue;
145  }
146  if (all_leafs[i_cnt]->data == NULL) {
147  some_thread_failed = -1;
148  continue;
149  }
150  if (((struct quaddata *)(all_leafs[i_cnt]->data))->points == NULL) {
151  continue;
152  }
153  else {
154  distx = (((struct quaddata *)(all_leafs[i_cnt]->data))->n_cols *
155  ew_res) *
156  0.1;
157  disty = (((struct quaddata *)(all_leafs[i_cnt]->data))->n_rows *
158  ns_res) *
159  0.1;
160  distxp = 0;
161  distyp = 0;
162  xmn = ((struct quaddata *)(all_leafs[i_cnt]->data))->x_orig;
163  xmx = ((struct quaddata *)(all_leafs[i_cnt]->data))->xmax;
164  ymn = ((struct quaddata *)(all_leafs[i_cnt]->data))->y_orig;
165  ymx = ((struct quaddata *)(all_leafs[i_cnt]->data))->ymax;
166  i = 0;
167  MAXENC = 0;
168  /* data is a window with zero points; some fields don't make
169  sense in this case so they are zero (like
170  resolution,dimensions */
171  /* CHANGE */
172  /* Calcutaing kmin for surrent segment (depends on the size) */
173 
174  /*****if (smseg <= 0.00001) MINPTS=params->kmin; else {} ***/
175  pr = pow(2., (xmx - xmn) / smseg - 1.);
176  MINPTS = params->kmin *
177  (pr / (1 + params->kmin * pr / params->KMAX2));
178  /* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf,
179  * smseg=%lf, DX=%lf \n",
180  * MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */
181 
182  data_local[tid] = (struct quaddata *)quad_data_new(
183  xmn - distx, ymn - disty, xmx + distx, ymx + disty, 0, 0, 0,
184  params->KMAX2);
185  npt = MT_region_data(info, tree, data_local[tid], params->KMAX2,
186  4);
187 
188  while ((npt < MINPTS) || (npt > params->KMAX2)) {
189  if (i >= 70) {
190  G_warning(_("Taking too long to find points for "
191  "interpolation - "
192  "please change the region to area where "
193  "your points are. "
194  "Continuing calculations..."));
195  break;
196  }
197  i++;
198  if (npt > params->KMAX2)
199  /* decrease window */
200  {
201  MAXENC = 1;
202  temp1 = distxp;
203  distxp = distx;
204  distx = distxp - fabs(distx - temp1) * 0.5;
205  temp2 = distyp;
206  distyp = disty;
207  disty = distyp - fabs(disty - temp2) * 0.5;
208  /* decrease by 50% of a previous change in window */
209  }
210  else {
211  temp1 = distyp;
212  distyp = disty;
213  temp2 = distxp;
214  distxp = distx;
215  if (MAXENC) {
216  disty = fabs(disty - temp1) * 0.5 + distyp;
217  distx = fabs(distx - temp2) * 0.5 + distxp;
218  }
219  else {
220  distx += distx;
221  disty += disty;
222  }
223  /* decrease by 50% of extra distance */
224  }
225  data_local[tid]->x_orig = xmn - distx; /* update window */
226  data_local[tid]->y_orig = ymn - disty;
227  data_local[tid]->xmax = xmx + distx;
228  data_local[tid]->ymax = ymx + disty;
229  data_local[tid]->n_points = 0;
230  npt = MT_region_data(info, tree, data_local[tid],
231  params->KMAX2, 4);
232  }
233 
234  if (totsegm != 0 && tid == 0) {
235  G_percent(cursegm, totsegm, 1);
236  }
237  data_local[tid]->n_rows =
238  ((struct quaddata *)(all_leafs[i_cnt]->data))->n_rows;
239  data_local[tid]->n_cols =
240  ((struct quaddata *)(all_leafs[i_cnt]->data))->n_cols;
241 
242  /* for printing out overlapping segments */
243  ((struct quaddata *)(all_leafs[i_cnt]->data))->x_orig =
244  xmn - distx;
245  ((struct quaddata *)(all_leafs[i_cnt]->data))->y_orig =
246  ymn - disty;
247  ((struct quaddata *)(all_leafs[i_cnt]->data))->xmax =
248  xmx + distx;
249  ((struct quaddata *)(all_leafs[i_cnt]->data))->ymax =
250  ymx + disty;
251 
252  data_local[tid]->x_orig = xmn;
253  data_local[tid]->y_orig = ymn;
254  data_local[tid]->xmax = xmx;
255  data_local[tid]->ymax = ymx;
256 
257  /* allocate memory for CV points only if cv is performed */
258  if (params->cv) {
259  if (!(point = (struct triple *)G_malloc(
260  sizeof(struct triple) *
261  data_local[tid]->n_points))) {
262  G_warning(_("Out of memory"));
263  some_thread_failed = -1;
264  continue;
265  }
266  }
267 
268  /*normalize the data so that the side of average segment is
269  * about 1m */
270  /* put data_points into point only if CV is performed */
271 
272  for (i = 0; i < data_local[tid]->n_points; i++) {
273  data_local[tid]->points[i].x =
274  (data_local[tid]->points[i].x -
275  data_local[tid]->x_orig) /
276  dnorm;
277  data_local[tid]->points[i].y =
278  (data_local[tid]->points[i].y -
279  data_local[tid]->y_orig) /
280  dnorm;
281  if (params->cv) {
282  point[i].x = data_local[tid]->points[i].x; /*cv stuff */
283  point[i].y = data_local[tid]->points[i].y; /*cv stuff */
284  point[i].z = data_local[tid]->points[i].z; /*cv stuff */
285  }
286 
287  /* commented out by Helena january 1997 as this is not
288  necessary although it may be useful to put normalization
289  of z back? data->points[i].z = data->points[i].z / dnorm;
290  this made smoothing self-adjusting based on dnorm
291  if (params->rsm < 0.) data->points[i].sm =
292  data->points[i].sm / dnorm;
293  */
294  }
295 
296  target_point.x = 0;
297  target_point.y = 0;
298  target_point.z = 0;
299 
300  /* one time interpolation and devi */
301  if (!params->cv) {
302  if (/* params */
303  IL_matrix_create_alloc(params, data_local[tid]->points,
304  data_local[tid]->n_points,
305  matrix[tid], indx[tid],
306  A[tid]) < 0) {
307  some_thread_failed = -1;
308  continue;
309  }
310 
311  for (i = 0; i < data_local[tid]->n_points; i++) {
312  b[tid][i + 1] = data_local[tid]->points[i].z;
313  }
314  b[tid][0] = 0.;
315  G_lubksb(matrix[tid], data_local[tid]->n_points + 1,
316  indx[tid], b[tid]);
317  /* put here condition to skip error if not needed */
318 
319  if (!params->create_devi) { /* check_points only for one
320  time interpolation */
321  params->check_points(params, data_local[tid], b[tid],
322  ertot, zmin, dnorm, &target_point);
323  }
324  }
325 
326  npoints = (params->cv || params->create_devi)
327  ? data_local[tid]->n_points
328  : 0;
329  for (point_index = 0; point_index < npoints;
330  point_index++) { /* loop only for cv or devi*/
331  if (params->cv) { /* cv: skip one point */
332  /* skip point for cv */
333  target_point.x = point[point_index].x;
334  target_point.y = point[point_index].y;
335  target_point.z = point[point_index].z;
336 
337  xx = target_point.x * dnorm + data_local[tid]->x_orig +
338  params->x_orig;
339  yy = target_point.y * dnorm + data_local[tid]->y_orig +
340  params->y_orig;
341  /* zz = point[point_index].z; */
342 
343  if (xx >= data_local[tid]->x_orig + params->x_orig &&
344  xx <= data_local[tid]->xmax + params->x_orig &&
345  yy >= data_local[tid]->y_orig + params->y_orig &&
346  yy <= data_local[tid]->ymax + params->y_orig) {
347 
348  j = 0;
349  for (k = 0; k < npoints; k++) {
350  if (k != point_index) {
351  data_local[tid]->points[j].x = point[k].x;
352  data_local[tid]->points[j].y = point[k].y;
353  data_local[tid]->points[j].z = point[k].z;
354  j++;
355  }
356  }
357 
358  /* segment area test for cv */
359  if (/* params */
361  params, data_local[tid]->points,
362  data_local[tid]->n_points - 1, matrix[tid],
363  indx[tid], A[tid]) < 0) {
364  some_thread_failed = -1;
365  continue;
366  }
367 
368  for (i = 0; i < data_local[tid]->n_points - 1;
369  i++) {
370  b[tid][i + 1] = data_local[tid]->points[i].z;
371  }
372  b[tid][0] = 0.;
373  G_lubksb(matrix[tid], data_local[tid]->n_points,
374  indx[tid], b[tid]);
375  }
376  else {
377  continue;
378  }
379  } /* cv: skip one point */
380 
381  if (params->create_devi) {
382  target_point.x = data_local[tid]->points[point_index].x;
383  target_point.y = data_local[tid]->points[point_index].y;
384  target_point.z = data_local[tid]->points[point_index].z;
385  }
386 
387  /* x, y, z is required input, while xmm, ymm, err output*/
388  pointz = target_point.z;
389  params->check_points(params, data_local[tid], b[tid], ertot,
390  zmin, dnorm, &target_point);
391 
392  err = target_point.z;
393  target_point.z = pointz + zmin;
394  xmm = target_point.x;
395  ymm = target_point.y;
396 
397  /* write out vector (point), if the point is inside the
398  * region*/
399  if (xmm >= data_local[tid]->x_orig + params->x_orig &&
400  xmm <= data_local[tid]->xmax + params->x_orig &&
401  ymm >= data_local[tid]->y_orig + params->y_orig &&
402  ymm <= data_local[tid]->ymax + params->y_orig) {
403  /* vect append, count, vect_write, db_execute will have
404  * conflicts between threads */
405 #pragma omp critical
406  {
407  params->check_points(params, NULL, NULL, &err, 0.0,
408  0.0, &target_point);
409  }
410  }
411  } /* end of computations for every point in cv or devi*/
412 
413  /* write out grid*/
414  if (!params->cv) {
415  if ((params->Tmp_fd_z != NULL) ||
416  (params->Tmp_fd_dx != NULL) ||
417  (params->Tmp_fd_dy != NULL) ||
418  (params->Tmp_fd_xx != NULL) ||
419  (params->Tmp_fd_yy != NULL) ||
420  (params->Tmp_fd_xy != NULL)) {
421 #pragma omp critical
422  {
423  if (params->grid_calc(params, data_local[tid],
424  bitmask, zmin, zmax, zminac,
425  zmaxac, gmin, gmax, c1min,
426  c1max, c2min, c2max, ertot,
427  b[tid], offset1, dnorm) < 0) {
428  some_thread_failed = -1;
429  }
430  }
431  }
432  }
433 
434  /* show after to catch 100% */
435 #pragma omp atomic
436  cursegm++;
437  if (totsegm < cursegm) {
438  G_debug(1, "%d %d", totsegm, cursegm);
439  }
440 
441  if (totsegm != 0 && tid == 0) {
442  G_percent(cursegm, totsegm, 1);
443  }
444  /*
445  G_free_matrix(matrix);
446  G_free_ivector(indx);
447  G_free_vector(b);
448  */
449  G_free(data_local[tid]->points);
450  G_free(data_local[tid]);
451  }
452  }
453  } /* All threads join master thread and terminate */
454 
455  for (i_cnt = 0; i_cnt < threads; i_cnt++) {
456  G_free(matrix[i_cnt]);
457  G_free(indx[i_cnt]);
458  G_free(b[i_cnt]);
459  G_free(A[i_cnt]);
460  }
461  G_free(all_leafs);
462  G_free(data_local);
463  G_free(matrix);
464  G_free(indx);
465  G_free(b);
466  G_free(A);
467 
468  if (some_thread_failed != 0) {
469  return -1;
470  }
471  return 1;
472 }
473 
474 /* cut given tree into separate leafs */
475 int cut_tree(struct multtree *tree, /* tree we want to cut */
476  struct multtree **cut_leafs, /* array of leafs */
477  int *where_to_add /* index of leaf which will be next */)
478 {
479  if (tree == NULL)
480  return -1;
481  if (tree->data == NULL)
482  return -1;
483  if (((struct quaddata *)(tree->data))->points == NULL) {
484  int i;
485 
486  for (i = 0; i < 4; i++) {
487  cut_tree(tree->leafs[i], cut_leafs, where_to_add);
488  }
489  return 1;
490  }
491  else {
492  cut_leafs[*where_to_add] = tree;
493  (*where_to_add)++;
494  return 1;
495  }
496 }
#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 void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
#define G_malloc(n)
Definition: defs/gis.h:94
void G_message(const char *,...) __attribute__((format(printf
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 IL_matrix_create_alloc(struct interp_params *, struct triple *, int, double **, int *, double *)
Creates system of linear equations from interpolated points.
Definition: matrix.c:69
double smallest_segment(struct multtree *, int)
Definition: segmen2d.c:338
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
int IL_interp_segments_2d_parallel(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, int threads)
Definition: bitmap.h:17
check_points_fn * check_points
Definition: interpf.h:127
FILE * Tmp_fd_xx
Definition: interpf.h:118
FILE * Tmp_fd_xy
Definition: interpf.h:118
FILE * Tmp_fd_yy
Definition: interpf.h:118
grid_calc_fn * grid_calc
Definition: interpf.h:123
double x_orig
Definition: interpf.h:107
FILE * Tmp_fd_dx
Definition: interpf.h:118
double y_orig
Definition: interpf.h:107
FILE * Tmp_fd_z
Definition: interpf.h:118
FILE * Tmp_fd_dy
Definition: interpf.h:118
bool create_devi
Definition: interpf.h:121
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
double z
Definition: dataquad.h:41
double x
Definition: dataquad.h:39
double y
Definition: dataquad.h:40
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:216