21 #include <grass/gis.h>
22 #include <grass/glocale.h>
23 #include <grass/interpf.h>
24 #include <grass/gmath.h>
26 static double smallest_segment(
struct multtree *,
int);
53 double zmin,
double zmax,
54 double *zminac,
double *zmaxac,
55 double *gmin,
double *gmax,
56 double *c1min,
double *c1max,
57 double *c2min,
double *c2max,
63 double xmn, xmx, ymn, ymx, distx, disty, distxp, distyp, temp1, temp2;
64 int i, npt, nptprev, MAXENC;
66 static int cursegm = 0;
67 static double *
b =
NULL;
68 static int *indx =
NULL;
69 static double **matrix =
NULL;
70 double ew_res, ns_res;
71 static int first_time = 1;
77 int m_skip, skip_index, j, k, segtest;
82 smseg = smallest_segment(info->
root, 4);
97 for (i = 0; i < 4; i++) {
99 bitmask, zmin, zmax, zminac, zmaxac, gmin,
100 gmax, c1min, c1max, c2min, c2max, ertot,
101 totsegm, offset1, dnorm);
122 pr = pow(2., (xmx - xmn) / smseg - 1.);
124 params->
kmin * (pr / (1 + params->
kmin * pr / params->
KMAX2));
129 xmx + distx, ymx + disty, 0, 0,
133 while ((npt < MINPTS) || (npt > params->
KMAX2)) {
135 G_warning(_(
"Taking too long to find points for interpolation - "
136 "please change the region to area where your points are. "
137 "Continuing calculations..."));
141 if (npt > params->
KMAX2)
148 distx = distxp - fabs(distx - temp1) * 0.5;
151 disty = distyp - fabs(disty - temp2) * 0.5;
161 disty = fabs(disty - temp1) * 0.5 + distyp;
162 distx = fabs(distx - temp2) * 0.5 + distxp;
170 data->
x_orig = xmn - distx;
171 data->
y_orig = ymn - disty;
172 data->
xmax = xmx + distx;
173 data->
ymax = ymx + disty;
229 for (i = 0; i < data->
n_points; i++) {
259 for (skip_index = 0; skip_index < m_skip; skip_index++) {
263 xx = point[skip_index].
x * dnorm + data->
x_orig +
265 yy = point[skip_index].
y * dnorm + data->
y_orig +
267 zz = point[skip_index].
z;
273 skip_point.
x = point[skip_index].
x;
274 skip_point.
y = point[skip_index].
y;
275 skip_point.
z = point[skip_index].
z;
276 for (k = 0; k < m_skip; k++) {
277 if (k != skip_index && params->
cv) {
292 else if (segtest == 1) {
299 for (i = 0; i < data->
n_points; i++)
304 params->
check_points(params, data, b, ertot, zmin, dnorm,
307 else if (segtest == 1) {
308 for (i = 0; i < data->
n_points - 1; i++)
312 params->
check_points(params, data, b, ertot, zmin, dnorm,
322 if (params->
grid_calc(params, data, bitmask,
323 zmin, zmax, zminac, zmaxac, gmin, gmax,
324 c1min, c1max, c2min, c2max, ertot, b,
331 if (totsegm < cursegm)
332 G_debug(1,
"%d %d", totsegm, cursegm);
348 static double smallest_segment(
struct multtree *tree,
int n_leafs)
350 static int first_time = 1;
352 static double minside;
360 for (ii = 0; ii < n_leafs; ii++) {
361 side = smallest_segment(tree->
leafs[ii], n_leafs);
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
int IL_interp_segments_2d(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)
struct quaddata * quad_data_new(double x_or, double y_or, double xmax, double ymax, int rows, int cols, int n_points, int kmax)
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
check_points_fn * check_points
int G_debug(int level, const char *msg,...)
Print debugging message.
void G_percent(long n, long d, int s)
Print percent complete messages.
int MT_region_data(struct tree_info *info, struct multtree *tree, struct quaddata *data, int MAX, int n_leafs)
int * G_alloc_ivector(size_t n)
Vector matrix memory allocation.
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
void G_free(void *buf)
Free allocated memory.
void G_warning(const char *msg,...)
Print a warning message to stderr.