35 double zmin,
double zmax,
36 double *zminac,
double *zmaxac,
37 double *gmin,
double *gmax,
38 double *c1min,
double *c1max,
double *c2min,
42 double *dnorm,
int overlap,
int inp_rows,
int inp_cols,
int fdsmooth,
43 int fdinp,
double ns_res,
double ew_res,
double inp_ns_res,
44 double inp_ew_res,
int dtens)
47 int i, j, k,
l, m, m1, i1;
51 double x_or, y_or, xm, ym;
52 static int first = 1, new_first = 1;
56 int out_check_rows, out_check_cols;
57 int first_row, last_row;
58 int first_col, last_col;
61 int rem_out_row, rem_out_col;
62 int inp_seg_r, inp_seg_c,
81 xmax = xmin + ew_res * params->
nsizc;
82 ymax = ymin + ns_res * params->
nsizr;
83 prev = inp_rows * inp_cols;
84 if (prev <= params->kmax)
92 if (num < params->kmin) {
93 if (((params->
kmin - num) > (prev + 1 - params->
kmax)) &&
94 (prev + 1 < params->
KMAX2)) {
103 if ((num > params->
kmin) && (num + 1 < params->
kmax)) {
110 out_seg_r = params->
nsizr / div;
111 out_seg_c = params->
nsizc / div;
112 inp_seg_r = inp_rows / div;
113 inp_seg_c = inp_cols / div;
114 rem_out_col = params->
nsizc % div;
115 rem_out_row = params->
nsizr % div;
116 overlap1 =
min1(overlap, inp_seg_c - 1);
117 overlap1 =
min1(overlap1, inp_seg_r - 1);
122 p_size = inp_seg_c * inp_seg_r;
125 p_size = (overlap1 * 2 + inp_seg_c) * (overlap1 * 2 + inp_seg_r);
130 fprintf(stderr,
"Cannot allocate memory for in_points\n");
136 sqrt(((xmax - xmin) * (ymax - ymin) * p_size) / (inp_rows * inp_cols));
139 params->
fi = params->
fi * (*dnorm) / 1000.;
140 fprintf(stderr,
"dnorm = %f, rescaled tension = %f\n", *dnorm,
148 input_data(params, 1, inp_rows, in_points, fdsmooth, fdinp, inp_rows,
149 inp_cols, zmin, inp_ns_res, inp_ew_res);
153 xm = params->
nsizc * ew_res;
154 ym = params->
nsizr * ns_res;
159 for (k = 1; k <= p_size; k++) {
160 if (!Rast_is_f_null_value(&(in_points[k - 1].z))) {
161 data->points[m1].x = in_points[k - 1].
x / (*dnorm);
162 data->points[m1].y = in_points[k - 1].
y / (*dnorm);
165 data->points[m1].z = (double)(in_points[k - 1].z);
166 data->points[m1].sm = in_points[k - 1].
smooth;
173 fprintf(stderr,
"Cannot allocate memory for indx\n");
177 fprintf(stderr,
"Cannot allocate memory for matrix\n");
181 fprintf(stderr,
"Cannot allocate memory for b\n");
185 if (params->
matrix_create(params, data->points, m1, matrix, indx) < 0)
187 for (i = 0; i < m1; i++) {
188 b[i + 1] = data->points[i].z;
195 if (params->
grid_calc(params, data, bitmask, zmin, zmax, zminac, zmaxac,
196 gmin, gmax, c1min, c1max, c2min, c2max, ertot,
b,
197 offset1, *dnorm) < 0) {
198 fprintf(stderr,
"interpolation failed\n");
209 fprintf(stderr,
"dnorm in ressegm after grid before out= %f \n",
215 out_seg_r = params->
nsizr / div;
216 out_seg_c = params->
nsizc / div;
217 inp_seg_r = inp_rows / div;
218 inp_seg_c = inp_cols / div;
219 rem_out_col = params->
nsizc % div;
220 rem_out_row = params->
nsizr % div;
221 overlap1 =
min1(overlap, inp_seg_c - 1);
222 overlap1 =
min1(overlap1, inp_seg_r - 1);
229 for (i = 1; i <= div; i++) {
230 if (i <= div - rem_out_row)
236 ngstr = out_check_rows + 1;
237 nszr = ngstr +
n_rows - 1;
238 y_or = (ngstr - 1) * ns_res;
243 first_row = (int)(y_or / inp_ns_res) + 1;
244 if (first_row > overlap1) {
245 first_row -= overlap1;
246 last_row = first_row + inp_seg_r + overlap1 * 2 - 1;
247 if (last_row > inp_rows) {
248 first_row -= (last_row - inp_rows);
254 last_row = first_row + inp_seg_r + overlap1 * 2 - 1;
256 if ((last_row > inp_rows) || (first_row < 1)) {
257 fprintf(stderr,
"Row overlap too large!\n");
260 input_data(params, first_row, last_row, in_points, fdsmooth, fdinp,
261 inp_rows, inp_cols, zmin, inp_ns_res, inp_ew_res);
263 for (j = 1; j <= div; j++) {
264 if (j <= div - rem_out_col)
270 ngstc = out_check_cols + 1;
271 nszc = ngstc +
n_cols - 1;
272 x_or = (ngstc - 1) * ew_res;
274 first_col = (int)(x_or / inp_ew_res) + 1;
275 if (first_col > overlap1) {
276 first_col -= overlap1;
277 last_col = first_col + inp_seg_c + overlap1 * 2 - 1;
278 if (last_col > inp_cols) {
279 first_col -= (last_col - inp_cols);
285 last_col = first_col + inp_seg_c + overlap1 * 2 - 1;
287 if ((last_col > inp_cols) || (first_col < 1)) {
288 fprintf(stderr,
"Column overlap too large!\n");
297 x_or, y_or, xm, ym, nszr - ngstr + 1, nszc - ngstc + 1, 0,
301 for (k = 0; k <= last_row - first_row; k++) {
302 for (
l = first_col - 1;
l < last_col;
l++) {
303 index = k * inp_cols +
l;
304 if (!Rast_is_f_null_value(&(in_points[index].z))) {
307 if ((in_points[index].
x - x_or >= 0) &&
308 (in_points[index].y - y_or >= 0) &&
309 ((nszc - 1) * ew_res - in_points[index].
x >= 0) &&
310 ((nszr - 1) * ns_res - in_points[index].y >= 0))
313 (in_points[index].
x - x_or) / (*dnorm);
315 (in_points[index].
y - y_or) / (*dnorm);
318 data->points[m].z = (double)(in_points[index].z);
319 data->points[m].sm = in_points[index].
smooth;
331 if (m <= params->KMAX2)
334 data->n_points = params->
KMAX2;
336 cursegm = (i - 1) * div + j - 1;
347 write_zeros(params, data, offset1);
356 "Cannot allocate memory for b\n");
362 "Cannot allocate memory for new_indx\n");
366 params->
KMAX2 + 1))) {
368 "Cannot allocate memory for new_matrix\n");
373 data->n_points, new_matrix,
377 for (i1 = 0; i1 < m; i1++) {
378 b[i1 + 1] = data->points[i1].z;
381 G_lubksb(new_matrix, data->n_points + 1, new_indx,
b);
386 if (params->
grid_calc(params, data, bitmask, zmin, zmax,
387 zminac, zmaxac, gmin, gmax, c1min,
388 c1max, c2min, c2max, ertot,
b,
389 offset1, *dnorm) < 0) {
391 fprintf(stderr,
"interpolate() failed\n");
401 "Cannot allocate memory for b\n");
407 "Cannot allocate memory for indx\n");
411 params->
KMAX2 + 1))) {
413 "Cannot allocate memory for matrix\n");
418 data->n_points, matrix, indx) < 0)
421 for (i1 = 0; i1 < m; i1++)
422 b[i1 + 1] = data->points[i1].z;
424 G_lubksb(matrix, data->n_points + 1, indx,
b);
429 if (params->
grid_calc(params, data, bitmask, zmin, zmax,
430 zminac, zmaxac, gmin, gmax, c1min,
431 c1max, c2min, c2max, ertot,
b,
432 offset1, *dnorm) < 0) {
434 fprintf(stderr,
"interpolate() failed\n");
459 fprintf(stderr,
"dnorm in ressegm after grid before out2= %f \n", *dnorm);
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)