28 #include <grass/interpf.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;
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
int * G_alloc_ivector(size_t)
Vector matrix memory allocation.
double smallest_segment(struct multtree *, int)
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)
void G_free(void *)
Free allocated memory.
double ** G_alloc_matrix(int, int)
Matrix memory allocation.
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)
double * G_alloc_vector(size_t)
Vector matrix memory allocation.
check_points_fn * check_points
int IL_matrix_create_alloc(struct interp_params *, struct triple *, int, double **, int *, double *)
Creates system of linear equations from interpolated points.
void G_percent(long, long, int)
Print percent complete messages.
int MT_region_data(struct tree_info *info, struct multtree *tree, struct quaddata *data, int MAX, int n_leafs)
void G_warning(const char *,...) __attribute__((format(printf
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
int G_debug(int, const char *,...) __attribute__((format(printf