GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
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 
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <math.h>
23 #if defined(_OPENMP)
24 #include <omp.h>
25 #endif
26 #include <grass/gis.h>
27 #include <grass/glocale.h>
28 #include <grass/interpf.h>
29 #include <grass/gmath.h>
30 
31 static int cut_tree(struct multtree *, struct multtree **, int *);
32 
33 
34 /*!
35  * See documentation for IL_interp_segments_2d.
36  * This is a parallel processing implementation.
37  */
38 int IL_interp_segments_2d_parallel(struct interp_params *params, 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;
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 (!
76  (matrix[i_cnt] =
77  G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
78  G_fatal_error(_("Out of memory"));
79  return -1;
80  }
81  }
82 
83  for (i_cnt = 0; i_cnt < threads; i_cnt++) {
84  if (!(indx[i_cnt] = G_alloc_ivector(params->KMAX2 + 1))) {
85  G_fatal_error(_("Out of memory"));
86  return -1;
87  }
88  }
89 
90  for (i_cnt = 0; i_cnt < threads; i_cnt++) {
91  if (!(b[i_cnt] = G_alloc_vector(params->KMAX2 + 3))) {
92  G_fatal_error(_("Out of memory"));
93  return -1;
94  }
95  }
96 
97  for (i_cnt = 0; i_cnt < threads; i_cnt++) {
98  if (!
99  (A[i_cnt] =
100  G_alloc_vector((params->KMAX2 + 2) * (params->KMAX2 + 2) + 1))) {
101  G_fatal_error(_("Out of memory"));
102  return -1;
103  }
104  }
105 
106  smseg = smallest_segment(tree, 4);
107  cut_tree(tree, all_leafs, &i);
108 
109  G_message(_("Starting parallel work"));
110 #pragma omp parallel firstprivate(tid, i, j, zmin, zmax, tree, totsegm, offset1, dnorm, smseg, ertot, params, info, all_leafs, bitmask, b, indx, matrix, data_local, A) shared(cursegm, threads, some_thread_failed, zminac, zmaxac, gmin, gmax, c1min, c1max, c2min, c2max) default(none)
111  {
112 #pragma omp for schedule(dynamic)
113  for (i_cnt = 0; i_cnt < totsegm; i_cnt++) {
114  /* Obtain thread id */
115 #if defined(_OPENMP)
116  tid = omp_get_thread_num();
117 #endif
118 
119  double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1,
120  temp2;
121  int npt, nptprev, MAXENC;
122  double ew_res, ns_res;
123  int MINPTS;
124  double pr;
125  struct triple *point;
126  struct triple skip_point;
127  int m_skip, skip_index, k, segtest;
128  double xx, yy, zz;
129 
130 
131  //struct quaddata *data_local;
132 
133  ns_res = (((struct quaddata *)(tree->data))->ymax -
134  ((struct quaddata *)(tree->data))->y_orig) /
135  params->nsizr;
136  ew_res =
137  (((struct quaddata *)(tree->data))->xmax -
138  ((struct quaddata *)(tree->data))->x_orig) / params->nsizc;
139 
140  if (all_leafs[i_cnt] == NULL) {
141  some_thread_failed = -1;
142  continue;
143  }
144  if (all_leafs[i_cnt]->data == NULL) {
145  some_thread_failed = -1;
146  continue;
147  }
148  if (((struct quaddata *)(all_leafs[i_cnt]->data))->points == NULL) {
149  continue;
150  }
151  else {
152  distx =
153  (((struct quaddata *)(all_leafs[i_cnt]->data))->n_cols *
154  ew_res) * 0.1;
155  disty =
156  (((struct quaddata *)(all_leafs[i_cnt]->data))->n_rows *
157  ns_res) * 0.1;
158  distxp = 0;
159  distyp = 0;
160  xmn = ((struct quaddata *)(all_leafs[i_cnt]->data))->x_orig;
161  xmx = ((struct quaddata *)(all_leafs[i_cnt]->data))->xmax;
162  ymn = ((struct quaddata *)(all_leafs[i_cnt]->data))->y_orig;
163  ymx = ((struct quaddata *)(all_leafs[i_cnt]->data))->ymax;
164  i = 0;
165  MAXENC = 0;
166  /* data is a window with zero points; some fields don't make sense in this case
167  so they are zero (like resolution,dimensions */
168  /* CHANGE */
169  /* Calcutaing kmin for surrent segment (depends on the size) */
170 
171  /*****if (smseg <= 0.00001) MINPTS=params->kmin; else {} ***/
172  pr = pow(2., (xmx - xmn) / smseg - 1.);
173  MINPTS =
174  params->kmin * (pr /
175  (1 + params->kmin * pr / params->KMAX2));
176  /* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf, smseg=%lf, DX=%lf \n", MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */
177 
178  data_local[tid] =
179  (struct quaddata *)quad_data_new(xmn - distx, ymn - disty,
180  xmx + distx, ymx + disty,
181  0, 0, 0, params->KMAX2);
182  npt =
183  MT_region_data(info, tree, data_local[tid], params->KMAX2,
184  4);
185 
186  while ((npt < MINPTS) || (npt > params->KMAX2)) {
187  if (i >= 70) {
188  G_warning(_("Taking too long to find points for interpolation - "
189  "please change the region to area where your points are. "
190  "Continuing calculations..."));
191  break;
192  }
193  i++;
194  if (npt > params->KMAX2)
195  /* decrease window */
196  {
197  MAXENC = 1;
198  nptprev = npt;
199  temp1 = distxp;
200  distxp = distx;
201  distx = distxp - fabs(distx - temp1) * 0.5;
202  temp2 = distyp;
203  distyp = disty;
204  disty = distyp - fabs(disty - temp2) * 0.5;
205  /* decrease by 50% of a previous change in window */
206  }
207  else {
208  nptprev = npt;
209  temp1 = distyp;
210  distyp = disty;
211  temp2 = distxp;
212  distxp = distx;
213  if (MAXENC) {
214  disty = fabs(disty - temp1) * 0.5 + distyp;
215  distx = fabs(distx - temp2) * 0.5 + distxp;
216  }
217  else {
218  distx += distx;
219  disty += disty;
220  }
221  /* decrease by 50% of extra distance */
222  }
223  data_local[tid]->x_orig = xmn - distx; /* update window */
224  data_local[tid]->y_orig = ymn - disty;
225  data_local[tid]->xmax = xmx + distx;
226  data_local[tid]->ymax = ymx + disty;
227  data_local[tid]->n_points = 0;
228  npt =
229  MT_region_data(info, tree, data_local[tid],
230  params->KMAX2, 4);
231  }
232 
233  if (totsegm != 0 && tid == 0) {
234  G_percent(cursegm, totsegm, 1);
235  }
236  data_local[tid]->n_rows =
237  ((struct quaddata *)(all_leafs[i_cnt]->data))->n_rows;
238  data_local[tid]->n_cols =
239  ((struct quaddata *)(all_leafs[i_cnt]->data))->n_cols;
240 
241  /* for printing out overlapping segments */
242  ((struct quaddata *)(all_leafs[i_cnt]->data))->x_orig =
243  xmn - distx;
244  ((struct quaddata *)(all_leafs[i_cnt]->data))->y_orig =
245  ymn - disty;
246  ((struct quaddata *)(all_leafs[i_cnt]->data))->xmax =
247  xmx + distx;
248  ((struct quaddata *)(all_leafs[i_cnt]->data))->ymax =
249  ymx + disty;
250 
251  data_local[tid]->x_orig = xmn;
252  data_local[tid]->y_orig = ymn;
253  data_local[tid]->xmax = xmx;
254  data_local[tid]->ymax = ymx;
255 
256  /* allocate memory for CV points only if cv is performed */
257  if (params->cv) {
258  if (!
259  (point =
260  (struct triple *)G_malloc(sizeof(struct triple) *
261  data_local[tid]->
262  n_points))) {
263  G_warning(_("Out of memory"));
264  some_thread_failed = -1;
265  continue;
266  }
267  }
268 
269  /*normalize the data so that the side of average segment is 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) / dnorm;
276  data_local[tid]->points[i].y =
277  (data_local[tid]->points[i].y -
278  data_local[tid]->y_orig) / dnorm;
279  if (params->cv) {
280  point[i].x = data_local[tid]->points[i].x; /*cv stuff */
281  point[i].y = data_local[tid]->points[i].y; /*cv stuff */
282  point[i].z = data_local[tid]->points[i].z; /*cv stuff */
283  }
284 
285  /* commented out by Helena january 1997 as this is not necessary
286  although it may be useful to put normalization of z back?
287  data->points[i].z = data->points[i].z / dnorm;
288  this made smoothing self-adjusting based on dnorm
289  if (params->rsm < 0.) data->points[i].sm = data->points[i].sm / dnorm;
290  */
291  }
292 
293  /* cv stuff */
294  if (params->cv) {
295  m_skip = data_local[tid]->n_points;
296  }
297  else {
298  m_skip = 1;
299  }
300  /* remove after cleanup - this is just for testing */
301  skip_point.x = 0.;
302  skip_point.y = 0.;
303  skip_point.z = 0.;
304 
305  for (skip_index = 0; skip_index < m_skip; skip_index++) {
306  if (params->cv) {
307  segtest = 0;
308  j = 0;
309  xx = point[skip_index].x * dnorm +
310  data_local[tid]->x_orig + params->x_orig;
311  yy = point[skip_index].y * dnorm +
312  data_local[tid]->y_orig + params->y_orig;
313  zz = point[skip_index].z;
314  if (xx >= data_local[tid]->x_orig + params->x_orig &&
315  xx <= data_local[tid]->xmax + params->x_orig &&
316  yy >= data_local[tid]->y_orig + params->y_orig &&
317  yy <= data_local[tid]->ymax + params->y_orig) {
318  segtest = 1;
319  skip_point.x = point[skip_index].x;
320  skip_point.y = point[skip_index].y;
321  skip_point.z = point[skip_index].z;
322  for (k = 0; k < m_skip; k++) {
323  if (k != skip_index && params->cv) {
324  data_local[tid]->points[j].x = point[k].x;
325  data_local[tid]->points[j].y = point[k].y;
326  data_local[tid]->points[j].z = point[k].z;
327  j++;
328  }
329  }
330  } /* segment area test */
331  }
332  if (!params->cv) {
333  if ( /* params */
334  IL_matrix_create_alloc(params,
335  data_local[tid]->points,
336  data_local[tid]->
337  n_points, matrix[tid],
338  indx[tid],
339  A[tid]) < 0) {
340  some_thread_failed = -1;
341  continue;
342  }
343  }
344  else if (segtest == 1) {
345  if ( /* params */
346  IL_matrix_create_alloc(params,
347  data_local[tid]->points,
348  data_local[tid]->
349  n_points - 1,
350  matrix[tid], indx[tid],
351  A[tid]) < 0) {
352  some_thread_failed = -1;
353  continue;
354  }
355  }
356  if (!params->cv) {
357  for (i = 0; i < data_local[tid]->n_points; i++) {
358  b[tid][i + 1] = data_local[tid]->points[i].z;
359  }
360  b[tid][0] = 0.;
361  G_lubksb(matrix[tid], data_local[tid]->n_points + 1,
362  indx[tid], b[tid]);
363  /* put here condition to skip error if not needed */
364  params->check_points(params, data_local[tid], b[tid],
365  ertot, zmin, dnorm, skip_point);
366  }
367  else if (segtest == 1) {
368  for (i = 0; i < data_local[tid]->n_points - 1; i++) {
369  b[tid][i + 1] = data_local[tid]->points[i].z;
370  }
371  b[tid][0] = 0.;
372  G_lubksb(matrix[tid], data_local[tid]->n_points,
373  indx[tid], b[tid]);
374  params->check_points(params, data_local[tid], b[tid],
375  ertot, zmin, dnorm, skip_point);
376  }
377  } /*end of cv loop */
378 
379 
380  if (!params->cv) {
381  if ((params->Tmp_fd_z != NULL) ||
382  (params->Tmp_fd_dx != NULL) ||
383  (params->Tmp_fd_dy != NULL) ||
384  (params->Tmp_fd_xx != NULL) ||
385  (params->Tmp_fd_yy != NULL) ||
386  (params->Tmp_fd_xy != NULL)) {
387 #pragma omp critical
388  {
389  if (params->grid_calc
390  (params, data_local[tid], bitmask, zmin, zmax,
391  zminac, zmaxac, gmin, gmax, c1min, c1max,
392  c2min, c2max, ertot, b[tid], offset1,
393  dnorm) < 0) {
394  some_thread_failed = -1;
395  }
396  }
397  }
398  }
399 
400  /* show after to catch 100% */
401 #pragma omp atomic
402  cursegm++;
403  if (totsegm < cursegm) {
404  G_debug(1, "%d %d", totsegm, cursegm);
405  }
406 
407  if (totsegm != 0 && tid == 0) {
408  G_percent(cursegm, totsegm, 1);
409  }
410  /*
411  G_free_matrix(matrix);
412  G_free_ivector(indx);
413  G_free_vector(b);
414  */
415  G_free(data_local[tid]->points);
416  G_free(data_local[tid]);
417  }
418  }
419  } /* All threads join master thread and terminate */
420 
421  for (i_cnt = 0; i_cnt < threads; i_cnt++) {
422  G_free(matrix[i_cnt]);
423  G_free(indx[i_cnt]);
424  G_free(b[i_cnt]);
425  G_free(A[i_cnt]);
426  }
427  G_free(all_leafs);
428  G_free(data_local);
429  G_free(matrix);
430  G_free(indx);
431  G_free(b);
432  G_free(A);
433 
434  if (some_thread_failed != 0) {
435  return -1;
436  }
437  return 1;
438 }
439 
440 
441 /* cut given tree into separate leafs */
442 int cut_tree(struct multtree *tree, /* tree we want to cut */
443  struct multtree **cut_leafs, /* array of leafs */
444  int *where_to_add /* index of leaf which will be next */ )
445 {
446  if (tree == NULL)
447  return -1;
448  if (tree->data == NULL)
449  return -1;
450  if (((struct quaddata *)(tree->data))->points == NULL) {
451  int i;
452 
453  for (i = 0; i < 4; i++) {
454  cut_tree(tree->leafs[i], cut_leafs, where_to_add);
455  }
456  return 1;
457  }
458  else {
459  cut_leafs[*where_to_add] = tree;
460  (*where_to_add)++;
461  return 1;
462  }
463 }
#define G_malloc(n)
Definition: defs/gis.h:112
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
int * G_alloc_ivector(size_t)
Vector matrix memory allocation.
Definition: ialloc.c:41
double y_orig
Definition: dataquad.h:51
double smallest_segment(struct multtree *, int)
Definition: segmen2d.c:346
FILE * Tmp_fd_xy
Definition: interpf.h:92
FILE * Tmp_fd_yy
Definition: interpf.h:92
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)
double z
Definition: dataquad.h:44
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
Definition: bitmap.h:17
double y_orig
Definition: interpf.h:87
double ** G_alloc_matrix(int, int)
Matrix memory allocation.
Definition: dalloc.c:60
#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
void G_message(const char *,...) __attribute__((format(printf
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
grid_calc_fn * grid_calc
Definition: interpf.h:97
FILE * Tmp_fd_dx
Definition: interpf.h:92
double b
Definition: r_raster.c:39
double * G_alloc_vector(size_t)
Vector matrix memory allocation.
Definition: dalloc.c:41
check_points_fn * check_points
Definition: interpf.h:99
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:72
double x
Definition: dataquad.h:42
double x_orig
Definition: interpf.h:87
double ymax
Definition: dataquad.h:53
void G_percent(long, long, int)
Print percent complete messages.
Definition: percent.c:62
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
void G_warning(const char *,...) __attribute__((format(printf
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition: lu.c:104
#define _(str)
Definition: glocale.h:10
Definition: qtree.h:56
double xmax
Definition: dataquad.h:52
int n_rows
Definition: dataquad.h:54
int n_cols
Definition: dataquad.h:55
int G_debug(int, const char *,...) __attribute__((format(printf
double y
Definition: dataquad.h:43
struct multtree ** leafs
Definition: qtree.h:59
int n_points
Definition: dataquad.h:56