26 #include <grass/gis.h>
27 #include <grass/glocale.h>
28 #include <grass/interpf.h>
29 #include <grass/gmath.h>
41 double zmin,
double zmax,
42 double *zminac,
double *zmaxac,
43 double *gmin,
double *gmax,
44 double *c1min,
double *c1max,
45 double *c2min,
double *c2max,
49 double dnorm,
int threads)
51 int some_thread_failed = 0;
58 double ***matrix =
NULL;
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);
74 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
83 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
90 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
97 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
107 cut_tree(tree, all_leafs, &i);
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)
112 #pragma omp for schedule(dynamic)
113 for (i_cnt = 0; i_cnt < totsegm; i_cnt++) {
116 tid = omp_get_thread_num();
119 double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1,
121 int npt, nptprev, MAXENC;
122 double ew_res, ns_res;
127 int m_skip, skip_index, k, segtest;
140 if (all_leafs[i_cnt] ==
NULL) {
141 some_thread_failed = -1;
144 if (all_leafs[i_cnt]->data ==
NULL) {
145 some_thread_failed = -1;
161 xmx = ((
struct quaddata *)(all_leafs[i_cnt]->data))->
xmax;
163 ymx = ((
struct quaddata *)(all_leafs[i_cnt]->data))->
ymax;
172 pr = pow(2., (xmx - xmn) / smseg - 1.);
175 (1 + params->
kmin * pr / params->
KMAX2));
180 xmx + distx, ymx + disty,
181 0, 0, 0, params->
KMAX2);
186 while ((npt < MINPTS) || (npt > params->
KMAX2)) {
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..."));
194 if (npt > params->
KMAX2)
201 distx = distxp - fabs(distx - temp1) * 0.5;
204 disty = distyp - fabs(disty - temp2) * 0.5;
214 disty = fabs(disty - temp1) * 0.5 + distyp;
215 distx = fabs(distx - temp2) * 0.5 + distxp;
223 data_local[tid]->
x_orig = xmn - distx;
224 data_local[tid]->
y_orig = ymn - disty;
225 data_local[tid]->
xmax = xmx + distx;
226 data_local[tid]->
ymax = ymx + disty;
233 if (totsegm != 0 && tid == 0) {
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;
264 some_thread_failed = -1;
272 for (i = 0; i < data_local[tid]->
n_points; i++) {
274 (data_local[tid]->
points[i].
x -
275 data_local[tid]->
x_orig) / dnorm;
277 (data_local[tid]->
points[i].
y -
278 data_local[tid]->
y_orig) / dnorm;
280 point[i].
x = data_local[tid]->
points[i].
x;
281 point[i].
y = data_local[tid]->
points[i].
y;
282 point[i].
z = data_local[tid]->
points[i].
z;
305 for (skip_index = 0; skip_index < m_skip; skip_index++) {
309 xx = point[skip_index].
x * dnorm +
311 yy = point[skip_index].
y * dnorm +
313 zz = point[skip_index].
z;
315 xx <= data_local[tid]->
xmax + params->
x_orig &&
317 yy <= data_local[tid]->
ymax + params->
y_orig) {
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;
340 some_thread_failed = -1;
344 else if (segtest == 1) {
350 matrix[tid], indx[tid],
352 some_thread_failed = -1;
357 for (i = 0; i < data_local[tid]->
n_points; i++) {
358 b[tid][i + 1] = data_local[tid]->
points[i].
z;
365 ertot, zmin, dnorm, skip_point);
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;
375 ertot, zmin, dnorm, skip_point);
390 (params, data_local[tid], bitmask, zmin, zmax,
391 zminac, zmaxac, gmin, gmax, c1min, c1max,
392 c2min, c2max, ertot,
b[tid], offset1,
394 some_thread_failed = -1;
403 if (totsegm < cursegm) {
404 G_debug(1,
"%d %d", totsegm, cursegm);
407 if (totsegm != 0 && tid == 0) {
421 for (i_cnt = 0; i_cnt < threads; i_cnt++) {
434 if (some_thread_failed != 0) {
453 for (i = 0; i < 4; i++) {
454 cut_tree(tree->
leafs[i], cut_leafs, where_to_add);
459 cut_leafs[*where_to_add] = tree;