27 #include <grass/interpf.h> 32 double,
double,
double);
36 double zmin,
double zmax,
37 double *zminac,
double *zmaxac,
38 double *gmin,
double *gmax,
39 double *c1min,
double *c1max,
double *c2min,
double *c2max,
51 double inp_ew_res,
int dtens)
54 int i, j, k,
l, m, m1, i1;
57 int n_rows, n_cols, inp_r, inp_c;
58 double x_or, y_or, xm, ym;
59 static int first = 1, new_first = 1;
63 int inp_check_rows, inp_check_cols,
64 out_check_rows, out_check_cols;
65 int first_row, last_row;
66 int first_col, last_col;
69 int rem_out_row, rem_out_col;
70 int inp_seg_r, inp_seg_c,
89 xmax = xmin + ew_res * params->
nsizc;
90 ymax = ymin + ns_res * params->
nsizr;
91 prev = inp_rows * inp_cols;
92 if (prev <= params->kmax)
100 if (num < params->kmin) {
101 if (((params->
kmin - num) > (prev + 1 - params->
kmax)) &&
102 (prev + 1 < params->
KMAX2)) {
111 if ((num > params->
kmin) && (num + 1 < params->
kmax)) {
118 out_seg_r = params->
nsizr / div;
119 out_seg_c = params->
nsizc / div;
120 inp_seg_r = inp_rows / div;
121 inp_seg_c = inp_cols / div;
122 rem_out_col = params->
nsizc % div;
123 rem_out_row = params->
nsizr % div;
124 overlap1 =
min1(overlap, inp_seg_c - 1);
125 overlap1 =
min1(overlap1, inp_seg_r - 1);
132 p_size = inp_seg_c * inp_seg_r;
135 p_size = (overlap1 * 2 + inp_seg_c) * (overlap1 * 2 + inp_seg_r);
142 fprintf(stderr,
"Cannot allocate memory for in_points\n");
148 sqrt(((xmax - xmin) * (ymax -
149 ymin) * p_size) / (inp_rows * inp_cols));
152 params->
fi = params->
fi * (*dnorm) / 1000.;
153 fprintf(stderr,
"dnorm = %f, rescaled tension = %f\n", *dnorm,
161 input_data(params, 1, inp_rows, in_points, fdsmooth, fdinp, inp_rows,
162 inp_cols, zmin, inp_ns_res, inp_ew_res);
166 xm = params->
nsizc * ew_res;
167 ym = params->
nsizr * ns_res;
173 for (k = 1; k <= p_size; k++) {
175 data->
points[m1].
x = in_points[k - 1].
x / (*dnorm);
176 data->
points[m1].
y = in_points[k - 1].
y / (*dnorm);
178 data->
points[m1].
z = (double)(in_points[k - 1].z);
186 fprintf(stderr,
"Cannot allocate memory for indx\n");
190 fprintf(stderr,
"Cannot allocate memory for matrix\n");
194 fprintf(stderr,
"Cannot allocate memory for b\n");
200 for (i = 0; i < m1; i++) {
206 params->
check_points(params, data, b, ertot, zmin, *dnorm, triple);
208 if (params->
grid_calc(params, data, bitmask,
209 zmin, zmax, zminac, zmaxac, gmin, gmax,
210 c1min, c1max, c2min, c2max, ertot, b, offset1,
212 fprintf(stderr,
"interpolation failed\n");
223 fprintf(stderr,
"dnorm in ressegm after grid before out= %f \n",
229 out_seg_r = params->
nsizr / div;
230 out_seg_c = params->
nsizc / div;
231 inp_seg_r = inp_rows / div;
232 inp_seg_c = inp_cols / div;
233 rem_out_col = params->
nsizc % div;
234 rem_out_row = params->
nsizr % div;
235 overlap1 =
min1(overlap, inp_seg_c - 1);
236 overlap1 =
min1(overlap1, inp_seg_r - 1);
245 for (i = 1; i <= div; i++) {
246 if (i <= div - rem_out_row)
249 n_rows = out_seg_r + 1;
253 ngstr = out_check_rows + 1;
254 nszr = ngstr + n_rows - 1;
255 y_or = (ngstr - 1) * ns_res;
260 first_row = (int)(y_or / inp_ns_res) + 1;
261 if (first_row > overlap1) {
262 first_row -= overlap1;
263 last_row = first_row + inp_seg_r + overlap1 * 2 - 1;
264 if (last_row > inp_rows) {
265 first_row -= (last_row - inp_rows);
271 last_row = first_row + inp_seg_r + overlap1 * 2 - 1;
273 if ((last_row > inp_rows) || (first_row < 1)) {
274 fprintf(stderr,
"Row overlap too large!\n");
277 input_data(params, first_row, last_row, in_points, fdsmooth, fdinp,
278 inp_rows, inp_cols, zmin, inp_ns_res, inp_ew_res);
280 for (j = 1; j <= div; j++) {
281 if (j <= div - rem_out_col)
284 n_cols = out_seg_c + 1;
287 ngstc = out_check_cols + 1;
288 nszc = ngstc + n_cols - 1;
289 x_or = (ngstc - 1) * ew_res;
291 first_col = (int)(x_or / inp_ew_res) + 1;
292 if (first_col > overlap1) {
293 first_col -= overlap1;
294 last_col = first_col + inp_seg_c + overlap1 * 2 - 1;
295 if (last_col > inp_cols) {
296 first_col -= (last_col - inp_cols);
302 last_col = first_col + inp_seg_c + overlap1 * 2 - 1;
304 if ((last_col > inp_cols) || (first_col < 1)) {
305 fprintf(stderr,
"Column overlap too large!\n");
319 for (k = 0; k <= last_row - first_row; k++) {
320 for (l = first_col - 1; l < last_col; l++) {
321 index = k * inp_cols +
l;
324 if ((in_points[index].
x - x_or >= 0) &&
325 (in_points[index].
y - y_or >= 0) &&
326 ((nszc - 1) * ew_res - in_points[index].
x >= 0) &&
327 ((nszr - 1) * ns_res - in_points[index].
y >= 0))
330 (in_points[index].
x - x_or) / (*dnorm);
332 (in_points[index].
y - y_or) / (*dnorm);
334 data->
points[m].
z = (double)(in_points[index].z);
347 if (m <= params->KMAX2)
352 inp_check_cols += inp_c;
353 cursegm = (i - 1) * div + j - 1;
364 write_zeros(params, data, offset1);
373 "Cannot allocate memory for b\n");
379 "Cannot allocate memory for new_indx\n");
385 params->
KMAX2 + 1))) {
387 "Cannot allocate memory for new_matrix\n");
393 new_matrix, new_indx) < 0)
396 for (i1 = 0; i1 < m; i1++) {
397 b[i1 + 1] = data->
points[i1].
z;
405 if (params->
grid_calc(params, data, bitmask,
406 zmin, zmax, zminac, zmaxac, gmin,
407 gmax, c1min, c1max, c2min, c2max,
408 ertot, b, offset1, *dnorm) < 0) {
410 fprintf(stderr,
"interpolate() failed\n");
420 "Cannot allocate memory for b\n");
426 "Cannot allocate memory for indx\n");
432 params->
KMAX2 + 1))) {
434 "Cannot allocate memory for matrix\n");
443 for (i1 = 0; i1 < m; i1++)
444 b[i1 + 1] = data->
points[i1].
z;
451 if (params->
grid_calc(params, data, bitmask,
452 zmin, zmax, zminac, zmaxac, gmin,
453 gmax, c1min, c1max, c2min, c2max,
454 ertot, b, offset1, *dnorm) < 0) {
456 fprintf(stderr,
"interpolate() failed\n");
470 inp_check_rows += inp_r;
482 fprintf(stderr,
"dnorm in ressegm after grid before out2= %f \n", *dnorm);
489 int first_row,
int last_row,
491 int fdsmooth,
int fdinp,
492 int inp_rows,
int inp_cols,
493 double zmin,
double inp_ns_res,
double inp_ew_res)
506 for (m1 = 0; m1 <= last_row - first_row; m1++) {
511 y = params->
y_orig + (m1 + first_row - 1 + 0.5) * inp_ns_res;
512 for (m2 = 0; m2 < inp_cols; m2++) {
513 x = params->
x_orig + (m2 + 0.5) * inp_ew_res;
518 sm = (double)cellsmooth[m2];
522 points[m1 * inp_cols + m2].
x = x - params->
x_orig;
523 points[m1 * inp_cols + m2].
y = y - params->
y_orig;
525 points[m1 * inp_cols + m2].
z =
526 cellinp[m2] * params->
zmult - zmin;
534 points[m1 * inp_cols + m2].
smooth = sm;
550 double x_or = data->
x_orig;
551 double y_or = data->
y_orig;
556 int ngstc, nszc, ngstr, nszr;
557 off_t offset, offset2;
558 double ns_res, ew_res;
562 ew_res = (((
struct quaddata *)(data))->xmax -
567 cond1 = ((params->
adx !=
NULL) || (params->
ady !=
NULL) || cond2);
569 ngstc = (int)(x_or / ew_res + 0.5) + 1;
570 nszc = ngstc + n_cols - 1;
571 ngstr = (int)(y_or / ns_res + 0.5) + 1;
572 nszr = ngstr + n_rows - 1;
574 for (k = ngstr; k <= nszr; k++) {
575 offset = offset1 * (k - 1);
576 for (l = ngstc; l <= nszc; l++) {
598 offset2 = (offset + ngstc - 1) *
sizeof(
FCELL);
599 if (params->
wr_temp(params, ngstc, nszc, offset2) < 0)
int * G_alloc_ivector(size_t)
Vector matrix memory allocation.
int IL_resample_interp_segments_2d(struct interp_params *params, 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, off_t offset1, double *dnorm, int overlap, int inp_rows, int inp_cols, int fdsmooth, int fdinp, double ns_res, double ew_res, double inp_ns_res, double inp_ew_res, int dtens)
void G_free(void *)
Free allocated memory.
matrix_create_fn * matrix_create
double ** G_alloc_matrix(int, int)
Matrix memory allocation.
#define Rast_is_f_null_value(fcellVal)
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
void Rast_set_f_null_value(FCELL *, int)
To set a number of FCELL raster values to NULL.
FCELL * Rast_allocate_f_buf(void)
Allocates memory for a raster map of type FCELL.
void G_percent(long, long, int)
Print percent complete messages.
void G_lubksb(double **a, int n, int *indx, double b[])
LU backward substitution.
void Rast_set_d_null_value(DCELL *, int)
To set a number of DCELL raster values to NULL.
void Rast_get_f_row(int, FCELL *, int)
Get raster row (FCELL type)