Loading...
Searching...
No Matches
segmen2d_parallel.c
Go to the documentation of this file.
/*!
* \file segmen2d.c
*
* \author H. Mitasova, I. Kosinovsky, D. Gerdes (single core)
* \author Stanislav Zubal, Michal Lacko (OpenMP version)
* \author Anna Petrasova (OpenMP version GRASS integration)
*
* \copyright
* (C) 1993-2017 by Helena Mitasova and the GRASS Development Team
*
* \copyright
* This program is free software under the
* GNU General Public License (>=v2).
* Read the file COPYING that comes with GRASS
* for details.
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#if defined(_OPENMP)
#include <omp.h>
#endif
#include <grass/gis.h>
#include <grass/glocale.h>
#include <grass/interpf.h>
#include <grass/gmath.h>
static int cut_tree(struct multtree *, struct multtree **, int *);
/*!
* See documentation for IL_interp_segments_2d.
* This is a parallel processing implementation.
*/
int IL_interp_segments_2d_parallel(
struct interp_params *params,
struct tree_info *info, /*!< info for the quad tree */
struct multtree *tree, /*!< current leaf of the quad tree */
struct BM *bitmask, /*!< bitmask */
double zmin, double zmax, /*!< min and max input z-values */
double *zminac, double *zmaxac, /*!< min and max interp. z-values */
double *gmin, double *gmax, /*!< min and max inperp. slope val. */
double *c1min, double *c1max, /*!< min and max interp. curv. val. */
double *c2min, double *c2max, /*!< min and max interp. curv. val. */
double *ertot, /*!< total interplating func. error */
int totsegm, /*!< total number of segments */
off_t offset1, /*!< offset for temp file writing */
double dnorm, int threads)
{
int some_thread_failed = 0;
int tid = 0;
int i = 0;
int j = 0;
int i_cnt;
int cursegm = 0;
double smseg;
double ***matrix = NULL;
int **indx = NULL;
double **b = NULL;
double **A = NULL;
struct quaddata **data_local;
struct multtree **all_leafs;
all_leafs =
(struct multtree **)G_malloc(sizeof(struct multtree *) * totsegm);
data_local =
(struct quaddata **)G_malloc(sizeof(struct quaddata *) * threads);
matrix = (double ***)G_malloc(sizeof(double **) * threads);
indx = (int **)G_malloc(sizeof(int *) * threads);
b = (double **)G_malloc(sizeof(double *) * threads);
A = (double **)G_malloc(sizeof(double *) * threads);
for (i_cnt = 0; i_cnt < threads; i_cnt++) {
if (!(matrix[i_cnt] =
G_alloc_matrix(params->KMAX2 + 1, params->KMAX2 + 1))) {
G_fatal_error(_("Out of memory"));
return -1;
}
}
for (i_cnt = 0; i_cnt < threads; i_cnt++) {
if (!(indx[i_cnt] = G_alloc_ivector(params->KMAX2 + 1))) {
G_fatal_error(_("Out of memory"));
return -1;
}
}
for (i_cnt = 0; i_cnt < threads; i_cnt++) {
if (!(b[i_cnt] = G_alloc_vector(params->KMAX2 + 3))) {
G_fatal_error(_("Out of memory"));
return -1;
}
}
for (i_cnt = 0; i_cnt < threads; i_cnt++) {
if (!(A[i_cnt] = G_alloc_vector(
(params->KMAX2 + 2) * (params->KMAX2 + 2) + 1))) {
G_fatal_error(_("Out of memory"));
return -1;
}
}
smseg = smallest_segment(tree, 4);
cut_tree(tree, all_leafs, &i);
G_message(_("Starting parallel work"));
#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)
{
#pragma omp for schedule(dynamic)
for (i_cnt = 0; i_cnt < totsegm; i_cnt++) {
/* Obtain thread id */
#if defined(_OPENMP)
tid = omp_get_thread_num();
#endif
double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1,
temp2;
int npt, MAXENC;
double ew_res, ns_res;
int MINPTS;
double pr;
struct triple *point;
struct triple target_point;
int npoints, point_index, k;
double xx, yy /*, zz */;
double xmm, ymm, err, pointz;
// struct quaddata *data_local;
ns_res = (((struct quaddata *)(tree->data))->ymax -
((struct quaddata *)(tree->data))->y_orig) /
params->nsizr;
ew_res = (((struct quaddata *)(tree->data))->xmax -
((struct quaddata *)(tree->data))->x_orig) /
params->nsizc;
if (all_leafs[i_cnt] == NULL) {
some_thread_failed = -1;
continue;
}
if (all_leafs[i_cnt]->data == NULL) {
some_thread_failed = -1;
continue;
}
if (((struct quaddata *)(all_leafs[i_cnt]->data))->points == NULL) {
continue;
}
else {
distx = (((struct quaddata *)(all_leafs[i_cnt]->data))->n_cols *
ew_res) *
0.1;
disty = (((struct quaddata *)(all_leafs[i_cnt]->data))->n_rows *
ns_res) *
0.1;
distxp = 0;
distyp = 0;
xmn = ((struct quaddata *)(all_leafs[i_cnt]->data))->x_orig;
xmx = ((struct quaddata *)(all_leafs[i_cnt]->data))->xmax;
ymn = ((struct quaddata *)(all_leafs[i_cnt]->data))->y_orig;
ymx = ((struct quaddata *)(all_leafs[i_cnt]->data))->ymax;
i = 0;
MAXENC = 0;
/* data is a window with zero points; some fields don't make
sense in this case so they are zero (like
resolution,dimensions */
/* CHANGE */
/* Calcutaing kmin for surrent segment (depends on the size) */
/*****if (smseg <= 0.00001) MINPTS=params->kmin; else {} ***/
pr = pow(2., (xmx - xmn) / smseg - 1.);
MINPTS = params->kmin *
(pr / (1 + params->kmin * pr / params->KMAX2));
/* fprintf(stderr,"MINPTS=%d, KMIN=%d, KMAX=%d, pr=%lf,
* smseg=%lf, DX=%lf \n",
* MINPTS,params->kmin,params->KMAX2,pr,smseg,xmx-xmn); */
data_local[tid] = (struct quaddata *)quad_data_new(
xmn - distx, ymn - disty, xmx + distx, ymx + disty, 0, 0, 0,
params->KMAX2);
npt = MT_region_data(info, tree, data_local[tid], params->KMAX2,
4);
while ((npt < MINPTS) || (npt > params->KMAX2)) {
if (i >= 70) {
G_warning(_("Taking too long to find points for "
"interpolation - "
"please change the region to area where "
"your points are. "
"Continuing calculations..."));
break;
}
i++;
if (npt > params->KMAX2)
/* decrease window */
{
MAXENC = 1;
temp1 = distxp;
distxp = distx;
distx = distxp - fabs(distx - temp1) * 0.5;
temp2 = distyp;
distyp = disty;
disty = distyp - fabs(disty - temp2) * 0.5;
/* decrease by 50% of a previous change in window */
}
else {
temp1 = distyp;
distyp = disty;
temp2 = distxp;
distxp = distx;
if (MAXENC) {
disty = fabs(disty - temp1) * 0.5 + distyp;
distx = fabs(distx - temp2) * 0.5 + distxp;
}
else {
distx += distx;
disty += disty;
}
/* decrease by 50% of extra distance */
}
data_local[tid]->x_orig = xmn - distx; /* update window */
data_local[tid]->y_orig = ymn - disty;
data_local[tid]->xmax = xmx + distx;
data_local[tid]->ymax = ymx + disty;
data_local[tid]->n_points = 0;
npt = MT_region_data(info, tree, data_local[tid],
params->KMAX2, 4);
}
if (totsegm != 0 && tid == 0) {
G_percent(cursegm, totsegm, 1);
}
data_local[tid]->n_rows =
((struct quaddata *)(all_leafs[i_cnt]->data))->n_rows;
data_local[tid]->n_cols =
((struct quaddata *)(all_leafs[i_cnt]->data))->n_cols;
/* for printing out overlapping segments */
((struct quaddata *)(all_leafs[i_cnt]->data))->x_orig =
xmn - distx;
((struct quaddata *)(all_leafs[i_cnt]->data))->y_orig =
ymn - disty;
((struct quaddata *)(all_leafs[i_cnt]->data))->xmax =
xmx + distx;
((struct quaddata *)(all_leafs[i_cnt]->data))->ymax =
ymx + disty;
data_local[tid]->x_orig = xmn;
data_local[tid]->y_orig = ymn;
data_local[tid]->xmax = xmx;
data_local[tid]->ymax = ymx;
/* allocate memory for CV points only if cv is performed */
if (params->cv) {
if (!(point = (struct triple *)G_malloc(
sizeof(struct triple) *
data_local[tid]->n_points))) {
G_warning(_("Out of memory"));
some_thread_failed = -1;
continue;
}
}
/*normalize the data so that the side of average segment is
* about 1m */
/* put data_points into point only if CV is performed */
for (i = 0; i < data_local[tid]->n_points; i++) {
data_local[tid]->points[i].x =
(data_local[tid]->points[i].x -
data_local[tid]->x_orig) /
dnorm;
data_local[tid]->points[i].y =
(data_local[tid]->points[i].y -
data_local[tid]->y_orig) /
dnorm;
if (params->cv) {
point[i].x = data_local[tid]->points[i].x; /*cv stuff */
point[i].y = data_local[tid]->points[i].y; /*cv stuff */
point[i].z = data_local[tid]->points[i].z; /*cv stuff */
}
/* commented out by Helena january 1997 as this is not
necessary although it may be useful to put normalization
of z back? data->points[i].z = data->points[i].z / dnorm;
this made smoothing self-adjusting based on dnorm
if (params->rsm < 0.) data->points[i].sm =
data->points[i].sm / dnorm;
*/
}
target_point.x = 0;
target_point.y = 0;
target_point.z = 0;
/* one time interpolation and devi */
if (!params->cv) {
if (/* params */
IL_matrix_create_alloc(params, data_local[tid]->points,
data_local[tid]->n_points,
matrix[tid], indx[tid],
A[tid]) < 0) {
some_thread_failed = -1;
continue;
}
for (i = 0; i < data_local[tid]->n_points; i++) {
b[tid][i + 1] = data_local[tid]->points[i].z;
}
b[tid][0] = 0.;
G_lubksb(matrix[tid], data_local[tid]->n_points + 1,
indx[tid], b[tid]);
/* put here condition to skip error if not needed */
if (!params->create_devi) { /* check_points only for one
time interpolation */
params->check_points(params, data_local[tid], b[tid],
ertot, zmin, dnorm, &target_point);
}
}
npoints = (params->cv || params->create_devi)
? data_local[tid]->n_points
: 0;
for (point_index = 0; point_index < npoints;
point_index++) { /* loop only for cv or devi*/
if (params->cv) { /* cv: skip one point */
/* skip point for cv */
target_point.x = point[point_index].x;
target_point.y = point[point_index].y;
target_point.z = point[point_index].z;
xx = target_point.x * dnorm + data_local[tid]->x_orig +
params->x_orig;
yy = target_point.y * dnorm + data_local[tid]->y_orig +
params->y_orig;
/* zz = point[point_index].z; */
if (xx >= data_local[tid]->x_orig + params->x_orig &&
xx <= data_local[tid]->xmax + params->x_orig &&
yy >= data_local[tid]->y_orig + params->y_orig &&
yy <= data_local[tid]->ymax + params->y_orig) {
j = 0;
for (k = 0; k < npoints; k++) {
if (k != point_index) {
data_local[tid]->points[j].x = point[k].x;
data_local[tid]->points[j].y = point[k].y;
data_local[tid]->points[j].z = point[k].z;
j++;
}
}
/* segment area test for cv */
if (/* params */
IL_matrix_create_alloc(
params, data_local[tid]->points,
data_local[tid]->n_points - 1, matrix[tid],
indx[tid], A[tid]) < 0) {
some_thread_failed = -1;
continue;
}
for (i = 0; i < data_local[tid]->n_points - 1;
i++) {
b[tid][i + 1] = data_local[tid]->points[i].z;
}
b[tid][0] = 0.;
G_lubksb(matrix[tid], data_local[tid]->n_points,
indx[tid], b[tid]);
}
else {
continue;
}
} /* cv: skip one point */
if (params->create_devi) {
target_point.x = data_local[tid]->points[point_index].x;
target_point.y = data_local[tid]->points[point_index].y;
target_point.z = data_local[tid]->points[point_index].z;
}
/* x, y, z is required input, while xmm, ymm, err output*/
pointz = target_point.z;
params->check_points(params, data_local[tid], b[tid], ertot,
zmin, dnorm, &target_point);
err = target_point.z;
target_point.z = pointz + zmin;
xmm = target_point.x;
ymm = target_point.y;
/* write out vector (point), if the point is inside the
* region*/
if (xmm >= data_local[tid]->x_orig + params->x_orig &&
xmm <= data_local[tid]->xmax + params->x_orig &&
ymm >= data_local[tid]->y_orig + params->y_orig &&
ymm <= data_local[tid]->ymax + params->y_orig) {
/* vect append, count, vect_write, db_execute will have
* conflicts between threads */
#pragma omp critical
{
params->check_points(params, NULL, NULL, &err, 0.0,
0.0, &target_point);
}
}
} /* end of computations for every point in cv or devi*/
/* write out grid*/
if (!params->cv) {
if ((params->Tmp_fd_z != NULL) ||
(params->Tmp_fd_dx != NULL) ||
(params->Tmp_fd_dy != NULL) ||
(params->Tmp_fd_xx != NULL) ||
(params->Tmp_fd_yy != NULL) ||
(params->Tmp_fd_xy != NULL)) {
#pragma omp critical
{
if (params->grid_calc(params, data_local[tid],
bitmask, zmin, zmax, zminac,
zmaxac, gmin, gmax, c1min,
c1max, c2min, c2max, ertot,
b[tid], offset1, dnorm) < 0) {
some_thread_failed = -1;
}
}
}
}
/* show after to catch 100% */
#pragma omp atomic
cursegm++;
if (totsegm < cursegm) {
G_debug(1, "%d %d", totsegm, cursegm);
}
if (totsegm != 0 && tid == 0) {
G_percent(cursegm, totsegm, 1);
}
/*
G_free_matrix(matrix);
G_free_ivector(indx);
G_free_vector(b);
*/
G_free(data_local[tid]->points);
G_free(data_local[tid]);
}
}
} /* All threads join master thread and terminate */
for (i_cnt = 0; i_cnt < threads; i_cnt++) {
G_free(matrix[i_cnt]);
G_free(indx[i_cnt]);
G_free(b[i_cnt]);
G_free(A[i_cnt]);
}
G_free(all_leafs);
G_free(data_local);
G_free(matrix);
G_free(indx);
G_free(b);
G_free(A);
if (some_thread_failed != 0) {
return -1;
}
return 1;
}
/* cut given tree into separate leafs */
int cut_tree(struct multtree *tree, /* tree we want to cut */
struct multtree **cut_leafs, /* array of leafs */
int *where_to_add /* index of leaf which will be next */)
{
if (tree == NULL)
return -1;
if (tree->data == NULL)
return -1;
if (((struct quaddata *)(tree->data))->points == NULL) {
int i;
for (i = 0; i < 4; i++) {
cut_tree(tree->leafs[i], cut_leafs, where_to_add);
}
return 1;
}
else {
cut_leafs[*where_to_add] = tree;
(*where_to_add)++;
return 1;
}
}
G_free
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:147
NULL
#define NULL
Definition ccmath.h:32
G_alloc_matrix
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition dalloc.c:55
G_alloc_vector
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition dalloc.c:38
quad_data_new
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
G_debug
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
b
double b
Definition driver/set_window.c:5
G_fatal_error
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
G_message
void G_message(const char *msg,...)
Print a message to stderr.
Definition gis/error.c:89
G_warning
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
G_alloc_ivector
int * G_alloc_ivector(size_t n)
Vector matrix memory allocation.
Definition ialloc.c:38
IL_matrix_create_alloc
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
smallest_segment
double smallest_segment(struct multtree *, int)
Definition segmen2d.c:342
G_lubksb
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
Definition lu.c:103
G_percent
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition percent.c:61
MT_region_data
int MT_region_data(struct tree_info *info, struct multtree *tree, struct quaddata *data, int MAX, int n_leafs)
Definition qtree.c:186
IL_interp_segments_2d_parallel
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 segmen2d_parallel.c:36
interp_params
Definition interpf.h:70
interp_params::check_points
check_points_fn * check_points
Definition interpf.h:132
interp_params::Tmp_fd_xx
FILE * Tmp_fd_xx
Definition interpf.h:123
interp_params::Tmp_fd_xy
FILE * Tmp_fd_xy
Definition interpf.h:123
interp_params::Tmp_fd_yy
FILE * Tmp_fd_yy
Definition interpf.h:123
interp_params::grid_calc
grid_calc_fn * grid_calc
Definition interpf.h:128
interp_params::x_orig
double x_orig
Definition interpf.h:112
interp_params::cv
int cv
Definition interpf.h:114
interp_params::Tmp_fd_dx
FILE * Tmp_fd_dx
Definition interpf.h:123
interp_params::y_orig
double y_orig
Definition interpf.h:112
interp_params::nsizr
int nsizr
Definition interpf.h:92
interp_params::Tmp_fd_z
FILE * Tmp_fd_z
Definition interpf.h:123
interp_params::nsizc
int nsizc
Definition interpf.h:92
interp_params::Tmp_fd_dy
FILE * Tmp_fd_dy
Definition interpf.h:123
interp_params::kmin
int kmin
Definition interpf.h:84
interp_params::KMAX2
int KMAX2
Definition interpf.h:99
interp_params::create_devi
bool create_devi
Definition interpf.h:126
multtree
Definition qtree.h:52
multtree::leafs
struct multtree ** leafs
Definition qtree.h:54
multtree::data
struct quaddata * data
Definition qtree.h:53
quaddata
Definition dataquad.h:45
quaddata::ymax
double ymax
Definition dataquad.h:49
quaddata::y_orig
double y_orig
Definition dataquad.h:47
quaddata::x_orig
double x_orig
Definition dataquad.h:46
quaddata::points
struct triple * points
Definition dataquad.h:53
quaddata::n_points
int n_points
Definition dataquad.h:52
quaddata::xmax
double xmax
Definition dataquad.h:48
quaddata::n_cols
int n_cols
Definition dataquad.h:51
quaddata::n_rows
int n_rows
Definition dataquad.h:50
tree_info
Definition qtree.h:45
triple
Definition dataquad.h:38
triple::z
double z
Definition dataquad.h:41
triple::x
double x
Definition dataquad.h:39
triple::y
double y
Definition dataquad.h:40
err
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition symbol/read.c:216