GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
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
30static 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 skip_point;
129 int m_skip, skip_index, k, segtest;
130 double xx, yy /*, zz */;
131
132 // struct quaddata *data_local;
133
134 ns_res = (((struct quaddata *)(tree->data))->ymax -
135 ((struct quaddata *)(tree->data))->y_orig) /
136 params->nsizr;
137 ew_res = (((struct quaddata *)(tree->data))->xmax -
138 ((struct quaddata *)(tree->data))->x_orig) /
139 params->nsizc;
140
141 if (all_leafs[i_cnt] == NULL) {
142 some_thread_failed = -1;
143 continue;
144 }
145 if (all_leafs[i_cnt]->data == NULL) {
146 some_thread_failed = -1;
147 continue;
148 }
149 if (((struct quaddata *)(all_leafs[i_cnt]->data))->points == NULL) {
150 continue;
151 }
152 else {
153 distx = (((struct quaddata *)(all_leafs[i_cnt]->data))->n_cols *
154 ew_res) *
155 0.1;
156 disty = (((struct quaddata *)(all_leafs[i_cnt]->data))->n_rows *
157 ns_res) *
158 0.1;
159 distxp = 0;
160 distyp = 0;
161 xmn = ((struct quaddata *)(all_leafs[i_cnt]->data))->x_orig;
162 xmx = ((struct quaddata *)(all_leafs[i_cnt]->data))->xmax;
163 ymn = ((struct quaddata *)(all_leafs[i_cnt]->data))->y_orig;
164 ymx = ((struct quaddata *)(all_leafs[i_cnt]->data))->ymax;
165 i = 0;
166 MAXENC = 0;
167 /* data is a window with zero points; some fields don't make
168 sense in this case so they are zero (like
169 resolution,dimensions */
170 /* CHANGE */
171 /* Calcutaing kmin for surrent segment (depends on the size) */
172
173 /*****if (smseg <= 0.00001) MINPTS=params->kmin; else {} ***/
174 pr = pow(2., (xmx - xmn) / smseg - 1.);
175 MINPTS = params->kmin *
176 (pr / (1 + params->kmin * pr / params->KMAX2));
177 /* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf,
178 * smseg=%lf, DX=%lf \n",
179 * MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */
180
181 data_local[tid] = (struct quaddata *)quad_data_new(
182 xmn - distx, ymn - disty, xmx + distx, ymx + disty, 0, 0, 0,
183 params->KMAX2);
184 npt = MT_region_data(info, tree, data_local[tid], params->KMAX2,
185 4);
186
187 while ((npt < MINPTS) || (npt > params->KMAX2)) {
188 if (i >= 70) {
189 G_warning(_("Taking too long to find points for "
190 "interpolation - "
191 "please change the region to area where "
192 "your points are. "
193 "Continuing calculations..."));
194 break;
195 }
196 i++;
197 if (npt > params->KMAX2)
198 /* decrease window */
199 {
200 MAXENC = 1;
201 temp1 = distxp;
202 distxp = distx;
203 distx = distxp - fabs(distx - temp1) * 0.5;
204 temp2 = distyp;
205 distyp = disty;
206 disty = distyp - fabs(disty - temp2) * 0.5;
207 /* decrease by 50% of a previous change in window */
208 }
209 else {
210 temp1 = distyp;
211 distyp = disty;
212 temp2 = distxp;
213 distxp = distx;
214 if (MAXENC) {
215 disty = fabs(disty - temp1) * 0.5 + distyp;
216 distx = fabs(distx - temp2) * 0.5 + distxp;
217 }
218 else {
219 distx += distx;
220 disty += disty;
221 }
222 /* decrease by 50% of extra distance */
223 }
224 data_local[tid]->x_orig = xmn - distx; /* update window */
225 data_local[tid]->y_orig = ymn - disty;
226 data_local[tid]->xmax = xmx + distx;
227 data_local[tid]->ymax = ymx + disty;
228 data_local[tid]->n_points = 0;
229 npt = 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 (!(point = (struct triple *)G_malloc(
259 sizeof(struct triple) *
260 data_local[tid]->n_points))) {
261 G_warning(_("Out of memory"));
262 some_thread_failed = -1;
263 continue;
264 }
265 }
266
267 /*normalize the data so that the side of average segment is
268 * about 1m */
269 /* put data_points into point only if CV is performed */
270
271 for (i = 0; i < data_local[tid]->n_points; i++) {
272 data_local[tid]->points[i].x =
273 (data_local[tid]->points[i].x -
274 data_local[tid]->x_orig) /
275 dnorm;
276 data_local[tid]->points[i].y =
277 (data_local[tid]->points[i].y -
278 data_local[tid]->y_orig) /
279 dnorm;
280 if (params->cv) {
281 point[i].x = data_local[tid]->points[i].x; /*cv stuff */
282 point[i].y = data_local[tid]->points[i].y; /*cv stuff */
283 point[i].z = data_local[tid]->points[i].z; /*cv stuff */
284 }
285
286 /* commented out by Helena january 1997 as this is not
287 necessary although it may be useful to put normalization
288 of z back? data->points[i].z = data->points[i].z / dnorm;
289 this made smoothing self-adjusting based on dnorm
290 if (params->rsm < 0.) data->points[i].sm =
291 data->points[i].sm / dnorm;
292 */
293 }
294
295 /* cv stuff */
296 if (params->cv) {
297 m_skip = data_local[tid]->n_points;
298 }
299 else {
300 m_skip = 1;
301 }
302 /* remove after cleanup - this is just for testing */
303 skip_point.x = 0.;
304 skip_point.y = 0.;
305 skip_point.z = 0.;
306
307 for (skip_index = 0; skip_index < m_skip; skip_index++) {
308 if (params->cv) {
309 segtest = 0;
310 j = 0;
311 xx = point[skip_index].x * dnorm +
312 data_local[tid]->x_orig + params->x_orig;
313 yy = point[skip_index].y * dnorm +
314 data_local[tid]->y_orig + params->y_orig;
315 /* zz = point[skip_index].z; */
316 if (xx >= data_local[tid]->x_orig + params->x_orig &&
317 xx <= data_local[tid]->xmax + params->x_orig &&
318 yy >= data_local[tid]->y_orig + params->y_orig &&
319 yy <= data_local[tid]->ymax + params->y_orig) {
320 segtest = 1;
321 skip_point.x = point[skip_index].x;
322 skip_point.y = point[skip_index].y;
323 skip_point.z = point[skip_index].z;
324 for (k = 0; k < m_skip; k++) {
325 if (k != skip_index && params->cv) {
326 data_local[tid]->points[j].x = point[k].x;
327 data_local[tid]->points[j].y = point[k].y;
328 data_local[tid]->points[j].z = point[k].z;
329 j++;
330 }
331 }
332 } /* segment area test */
333 }
334 if (!params->cv) {
335 if (/* params */
337 params, data_local[tid]->points,
338 data_local[tid]->n_points, matrix[tid],
339 indx[tid], A[tid]) < 0) {
340 some_thread_failed = -1;
341 continue;
342 }
343 }
344 else if (segtest == 1) {
345 if (/* params */
347 params, data_local[tid]->points,
348 data_local[tid]->n_points - 1, matrix[tid],
349 indx[tid], A[tid]) < 0) {
350 some_thread_failed = -1;
351 continue;
352 }
353 }
354 if (!params->cv) {
355 for (i = 0; i < data_local[tid]->n_points; i++) {
356 b[tid][i + 1] = data_local[tid]->points[i].z;
357 }
358 b[tid][0] = 0.;
359 G_lubksb(matrix[tid], data_local[tid]->n_points + 1,
360 indx[tid], b[tid]);
361 /* put here condition to skip error if not needed */
362 params->check_points(params, data_local[tid], b[tid],
363 ertot, zmin, dnorm, skip_point);
364 }
365 else if (segtest == 1) {
366 for (i = 0; i < data_local[tid]->n_points - 1; i++) {
367 b[tid][i + 1] = data_local[tid]->points[i].z;
368 }
369 b[tid][0] = 0.;
370 G_lubksb(matrix[tid], data_local[tid]->n_points,
371 indx[tid], b[tid]);
372 params->check_points(params, data_local[tid], b[tid],
373 ertot, zmin, dnorm, skip_point);
374 }
375 } /*end of cv loop */
376
377 if (!params->cv) {
378 if ((params->Tmp_fd_z != NULL) ||
379 (params->Tmp_fd_dx != NULL) ||
380 (params->Tmp_fd_dy != NULL) ||
381 (params->Tmp_fd_xx != NULL) ||
382 (params->Tmp_fd_yy != NULL) ||
383 (params->Tmp_fd_xy != NULL)) {
384#pragma omp critical
385 {
386 if (params->grid_calc(params, data_local[tid],
387 bitmask, zmin, zmax, zminac,
388 zmaxac, gmin, gmax, c1min,
389 c1max, c2min, c2max, ertot,
390 b[tid], offset1, dnorm) < 0) {
391 some_thread_failed = -1;
392 }
393 }
394 }
395 }
396
397 /* show after to catch 100% */
398#pragma omp atomic
399 cursegm++;
400 if (totsegm < cursegm) {
401 G_debug(1, "%d %d", totsegm, cursegm);
402 }
403
404 if (totsegm != 0 && tid == 0) {
405 G_percent(cursegm, totsegm, 1);
406 }
407 /*
408 G_free_matrix(matrix);
409 G_free_ivector(indx);
410 G_free_vector(b);
411 */
412 G_free(data_local[tid]->points);
413 G_free(data_local[tid]);
414 }
415 }
416 } /* All threads join master thread and terminate */
417
418 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
419 G_free(matrix[i_cnt]);
420 G_free(indx[i_cnt]);
421 G_free(b[i_cnt]);
422 G_free(A[i_cnt]);
423 }
424 G_free(all_leafs);
425 G_free(data_local);
426 G_free(matrix);
427 G_free(indx);
428 G_free(b);
429 G_free(A);
430
431 if (some_thread_failed != 0) {
432 return -1;
433 }
434 return 1;
435}
436
437/* cut given tree into separate leafs */
438int cut_tree(struct multtree *tree, /* tree we want to cut */
439 struct multtree **cut_leafs, /* array of leafs */
440 int *where_to_add /* index of leaf which will be next */)
441{
442 if (tree == NULL)
443 return -1;
444 if (tree->data == NULL)
445 return -1;
446 if (((struct quaddata *)(tree->data))->points == NULL) {
447 int i;
448
449 for (i = 0; i < 4; i++) {
450 cut_tree(tree->leafs[i], cut_leafs, where_to_add);
451 }
452 return 1;
453 }
454 else {
455 cut_leafs[*where_to_add] = tree;
456 (*where_to_add)++;
457 return 1;
458 }
459}
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150
#define NULL
Definition ccmath.h:32
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition dalloc.c:57
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition dalloc.c:39
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
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
double b
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
void G_message(const char *msg,...)
Print a message to stderr.
Definition gis/error.c:89
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
int * G_alloc_ivector(size_t n)
Vector matrix memory allocation.
Definition ialloc.c:38
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
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition lu.c:103
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition percent.c:61
int MT_region_data(struct tree_info *info, struct multtree *tree, struct quaddata *data, int MAX, int n_leafs)
Definition qtree.c:186
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)
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
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