1
2 /**********************************************************************
3 *
4 * MODULE: r.resamp.bspline
5 *
6 * AUTHOR(S): Markus Metz
7 *
8 * PURPOSE: Spline Interpolation
9 *
10 * COPYRIGHT: (C) 2010, 2012 by GRASS development team
11 *
12 * This program is free software under the GNU General
13 * Public License (>=v2). Read the file COPYING that
14 * comes with GRASS for details.
15 *
16 **********************************************************************/
17
18 /* INCLUDES */
19 #include <stdlib.h>
20 #include <string.h>
21 #include <sys/types.h>
22 #include <sys/stat.h>
23 #include <fcntl.h>
24 #include <math.h>
25 #include "bspline.h"
26
27 #define SEGSIZE 64
28
29 /*--------------------------------------------------------------------*/
main(int argc,char * argv[])30 int main(int argc, char *argv[])
31 {
32 /* Variable declarations */
33 int nsply, nsplx, row, col, nrows, ncols, nsplx_adj, nsply_adj;
34 int nsubregion_col, nsubregion_row, subregion_row, subregion_col;
35 int subregion = 0, nsubregions = 0;
36 int last_row, last_column, interp_method; /* booleans */
37 double lambda, mean;
38 double N_extension, E_extension, edgeE, edgeN;
39 double stepN, stepE;
40
41 char title[64];
42
43 int dim_vect, nparameters, BW;
44 double *TN, *Q, *parVect; /* Interpolating and least-square vectors */
45 double **N, **obsVect; /* Interpolation and least-square matrix */
46
47 SEGMENT in_seg, out_seg, mask_seg;
48 char *in_file, *out_file, *mask_file;
49 double seg_size;
50 int seg_mb, segments_in_memory;
51 int have_mask;
52
53
54 int inrastfd, outrastfd;
55 DCELL *drastbuf, dval;
56 struct History history;
57
58 struct Map_info Grid;
59 struct line_pnts *Points;
60 struct line_cats *Cats;
61 int cat = 1;
62
63 struct GModule *module;
64 struct Option *in_opt, *out_opt, *grid_opt, *stepE_opt, *stepN_opt,
65 *lambda_f_opt, *method_opt, *mask_opt, *memory_opt;
66 struct Flag *null_flag, *cross_corr_flag;
67
68 struct Reg_dimens dims;
69 struct Cell_head elaboration_reg, src_reg, dest_reg;
70 struct bound_box general_box, overlap_box, dest_box;
71 struct bound_box last_overlap_box, last_general_box;
72
73 struct Point *observ;
74 struct Point *observ_marked;
75
76 /*----------------------------------------------------------------*/
77 /* Options declarations */
78 module = G_define_module();
79 G_add_keyword(_("raster"));
80 G_add_keyword(_("surface"));
81 G_add_keyword(_("resample"));
82 G_add_keyword(_("interpolation"));
83 G_add_keyword(_("splines"));
84 G_add_keyword(_("bilinear"));
85 G_add_keyword(_("bicubic"));
86 G_add_keyword(_("no-data filling"));
87 module->description =
88 _("Performs bilinear or bicubic spline interpolation with Tykhonov regularization.");
89
90 in_opt = G_define_standard_option(G_OPT_R_INPUT);
91
92 out_opt = G_define_standard_option(G_OPT_R_OUTPUT);
93
94 grid_opt = G_define_standard_option(G_OPT_V_OUTPUT);
95 grid_opt->key = "grid";
96 grid_opt->description = _("Name for output vector map with interpolation grid");
97 grid_opt->required = NO;
98
99 mask_opt = G_define_standard_option(G_OPT_R_INPUT);
100 mask_opt->key = "mask";
101 mask_opt->label = _("Name of raster map to use for masking");
102 mask_opt->description = _("Only cells that are not NULL and not zero are interpolated");
103 mask_opt->required = NO;
104
105 stepE_opt = G_define_option();
106 stepE_opt->key = "ew_step";
107 stepE_opt->type = TYPE_DOUBLE;
108 stepE_opt->required = NO;
109 stepE_opt->description =
110 _("Length of each spline step in the east-west direction. Default: 1.5 * ewres.");
111 stepE_opt->guisection = _("Settings");
112
113 stepN_opt = G_define_option();
114 stepN_opt->key = "ns_step";
115 stepN_opt->type = TYPE_DOUBLE;
116 stepN_opt->required = NO;
117 stepN_opt->description =
118 _("Length of each spline step in the north-south direction. Default: 1.5 * nsres.");
119 stepN_opt->guisection = _("Settings");
120
121 method_opt = G_define_standard_option(G_OPT_R_INTERP_TYPE);
122 method_opt->description = _("Spline interpolation algorithm");
123 method_opt->options = "bilinear,bicubic";
124 method_opt->answer = "bicubic";
125 method_opt->guisection = _("Settings");
126 G_asprintf((char **) &(method_opt->descriptions),
127 "bilinear;%s;bicubic;%s",
128 _("Bilinear interpolation"),
129 _("Bicubic interpolation"));
130
131 lambda_f_opt = G_define_option();
132 lambda_f_opt->key = "lambda";
133 lambda_f_opt->type = TYPE_DOUBLE;
134 lambda_f_opt->required = NO;
135 lambda_f_opt->description = _("Tykhonov regularization parameter (affects smoothing)");
136 lambda_f_opt->answer = "0.01";
137 lambda_f_opt->guisection = _("Settings");
138
139 null_flag = G_define_flag();
140 null_flag->key = 'n';
141 null_flag->label = _("Only interpolate null cells in input raster map");
142 null_flag->guisection = _("Settings");
143
144 cross_corr_flag = G_define_flag();
145 cross_corr_flag->key = 'c';
146 cross_corr_flag->description =
147 _("Find the best Tykhonov regularizing parameter using a \"leave-one-out\" cross validation method");
148
149 memory_opt = G_define_standard_option(G_OPT_MEMORYMB);
150
151 /*----------------------------------------------------------------*/
152 /* Parsing */
153 G_gisinit(argv[0]);
154 if (G_parser(argc, argv))
155 exit(EXIT_FAILURE);
156
157 if (!strcmp(method_opt->answer, "bilinear"))
158 interp_method = P_BILINEAR;
159 else
160 interp_method = P_BICUBIC;
161
162 lambda = atof(lambda_f_opt->answer);
163
164 /* Setting regions and boxes */
165 G_debug(1, "Interpolation: Setting regions and boxes");
166 G_get_set_window(&dest_reg);
167 G_get_set_window(&elaboration_reg);
168 Vect_region_box(&dest_reg, &dest_box);
169 Vect_region_box(&elaboration_reg, &overlap_box);
170 Vect_region_box(&elaboration_reg, &general_box);
171
172 nrows = Rast_window_rows();
173 ncols = Rast_window_cols();
174
175 /* get window of input map */
176 Rast_get_cellhd(in_opt->answer, "", &src_reg);
177
178 if (stepE_opt->answer) {
179 stepE = atof(stepE_opt->answer);
180 if (stepE <= .0)
181 G_fatal_error(_("ew_step must be positive"));
182 }
183 else
184 stepE = src_reg.ew_res * 1.5;
185
186 if (stepN_opt->answer) {
187 stepN = atof(stepN_opt->answer);
188 if (stepN <= .0)
189 G_fatal_error(_("ns_step must be positive"));
190 }
191 else
192 stepN = src_reg.ns_res * 1.5;
193
194 /*------------------------------------------------------------------
195 | Subdividing and working with tiles:
196 | Each original region will be divided into several subregions.
197 | Each one will be overlaped by its neighbouring subregions.
198 | The overlapping is calculated as a fixed OVERLAP_SIZE times
199 | the largest spline step plus 2 * orlo
200 ----------------------------------------------------------------*/
201
202 /* Fixing parameters of the elaboration region */
203 P_zero_dim(&dims); /* Set dim struct to zero */
204
205 nsplx_adj = NSPLX_MAX;
206 nsply_adj = NSPLY_MAX;
207 if (interp_method == P_BICUBIC) {
208 nsplx_adj = 100;
209 nsply_adj = 100;
210 }
211
212 dims.overlap = OVERLAP_SIZE * (stepN > stepE ? stepN : stepE);
213 P_get_edge(interp_method, &dims, stepE, stepN);
214 P_set_dim(&dims, stepE, stepN, &nsplx_adj, &nsply_adj);
215
216 G_verbose_message(_("spline step in ew direction %g"), stepE);
217 G_verbose_message(_("spline step in ns direction %g"), stepN);
218 G_verbose_message(_("adjusted EW splines %d"), nsplx_adj);
219 G_verbose_message(_("adjusted NS splines %d"), nsply_adj);
220
221 /* calculate number of subregions */
222 edgeE = dims.ew_size - dims.overlap - 2 * dims.edge_v;
223 edgeN = dims.sn_size - dims.overlap - 2 * dims.edge_h;
224
225 N_extension = dest_reg.north - dest_reg.south;
226 E_extension = dest_reg.east - dest_reg.west;
227
228 nsubregion_col = ceil(E_extension / edgeE) + 0.5;
229 nsubregion_row = ceil(N_extension / edgeN) + 0.5;
230
231 if (nsubregion_col < 0)
232 nsubregion_col = 0;
233 if (nsubregion_row < 0)
234 nsubregion_row = 0;
235
236 nsubregions = nsubregion_row * nsubregion_col;
237
238 if (G_projection() == PROJECTION_LL) {
239 /* try to shift source window to overlap with destination window */
240 while (src_reg.west >= dest_reg.east && src_reg.east - 360.0 > dest_reg.west) {
241 src_reg.east -= 360.0;
242 src_reg.west -= 360.0;
243 }
244 while (src_reg.east <= dest_reg.west && src_reg.west + 360.0 < dest_reg.east) {
245 src_reg.east += 360.0;
246 src_reg.west += 360.0;
247 }
248 }
249
250 G_debug(1, "-------------------------------------");
251 G_debug(1, "source north %f", src_reg.north);
252 G_debug(1, "source south %f", src_reg.south);
253 G_debug(1, "source west %f", src_reg.west);
254 G_debug(1, "source east %f", src_reg.east);
255 G_debug(1, "-------------------------------------");
256
257 /* adjust source window */
258 {
259 double north = dest_reg.north + 2 * dims.edge_h;
260 double south = dest_reg.south - 2 * dims.edge_h;
261 int r0 = (int)(floor(Rast_northing_to_row(north, &src_reg)) - 0.5);
262 int r1 = (int)(ceil(Rast_northing_to_row(south, &src_reg)) + 0.5);
263 double east = dest_reg.east + 2 * dims.edge_v;
264 double west = dest_reg.west - 2 * dims.edge_v;
265 /* do not use Rast_easting_to_col() because of G_adjust_easting() */
266 /*
267 int c0 = (int)ceil(Rast_easting_to_col(east, &src_reg) + 0.5);
268 int c1 = (int)floor(Rast_easting_to_col(west, &src_reg) + 0.5);
269 */
270 int c0 = (int)(ceil(((east - src_reg.west) / src_reg.ew_res)) + 0.5);
271 int c1 = (int)(floor(((west - src_reg.west) / src_reg.ew_res)) - 0.5);
272
273 src_reg.north -= src_reg.ns_res * (r0);
274 src_reg.south -= src_reg.ns_res * (r1 - src_reg.rows);
275 src_reg.east += src_reg.ew_res * (c0 - src_reg.cols);
276 src_reg.west += src_reg.ew_res * (c1);
277 src_reg.rows = r1 - r0;
278 src_reg.cols = c0 - c1;
279 }
280
281 /* switch to buffered input raster window */
282 G_set_window(&src_reg);
283 Rast_set_window(&src_reg);
284
285 G_debug(1, "new source north %f", src_reg.north);
286 G_debug(1, "new source south %f", src_reg.south);
287 G_debug(1, "new source west %f", src_reg.west);
288 G_debug(1, "new source east %f", src_reg.east);
289 G_debug(1, "-------------------------------------");
290
291 nrows = Rast_window_rows();
292 ncols = Rast_window_cols();
293
294 G_debug(1, "%d new rows, %d new cols", nrows, ncols);
295
296 seg_mb = atoi(memory_opt->answer);
297 if (seg_mb < 3)
298 G_fatal_error(_("Memory in MB must be >= 3"));
299
300 if (mask_opt->answer || null_flag->answer) {
301 seg_size = sizeof(double) * 2 + sizeof(char);
302 }
303 else {
304 seg_size = sizeof(double) * 2;
305 }
306 seg_size = (seg_size * SEGSIZE * SEGSIZE) / (1 << 20);
307 segments_in_memory = seg_mb / seg_size + 0.5;
308 if (segments_in_memory < 1)
309 segments_in_memory = 1;
310
311 in_file = G_tempfile();
312 if (Segment_open(&in_seg, in_file, nrows, ncols, SEGSIZE, SEGSIZE,
313 sizeof(double), segments_in_memory) != 1)
314 G_fatal_error(_("Can not create temporary file"));
315
316 /* read raster input */
317 inrastfd = Rast_open_old(in_opt->answer, "");
318 drastbuf = Rast_allocate_buf(DCELL_TYPE);
319
320 G_message(_("Loading input raster <%s>"), in_opt->answer);
321 if (1) {
322 int got_one = 0;
323 for (row = 0; row < nrows; row++) {
324 G_percent(row, nrows, 9);
325
326 Rast_get_d_row_nomask(inrastfd, drastbuf, row);
327
328 for (col = 0; col < ncols; col++) {
329 dval = drastbuf[col];
330 if (!Rast_is_d_null_value(&dval)) {
331 got_one++;
332 }
333 if (Segment_put(&in_seg, &dval, row, col) < 0)
334 G_fatal_error(_("Can not write to temp file"));;
335 }
336
337 }
338 if (!got_one)
339 G_fatal_error(_("Only NULL cells in input raster"));
340 }
341 G_percent(row, nrows, 2);
342 G_free(drastbuf);
343 Rast_close(inrastfd);
344
345 /* switch back to destination = current window */
346 G_set_window(&dest_reg);
347 Rast_set_window(&dest_reg);
348 nrows = Rast_window_rows();
349 ncols = Rast_window_cols();
350
351 G_debug(1, "-------------------------------------");
352 G_debug(1, "dest north %f", dest_reg.north);
353 G_debug(1, "dest south %f", dest_reg.south);
354 G_debug(1, "dest west %f", dest_reg.west);
355 G_debug(1, "dest east %f", dest_reg.east);
356 G_debug(1, "-------------------------------------");
357
358 /* cross-correlation */
359 if (cross_corr_flag->answer) {
360 G_debug(1, "CrossCorrelation()");
361
362 if (cross_correlation(&in_seg, &src_reg, stepE, stepN) != TRUE)
363 G_fatal_error(_("Cross validation didn't finish correctly"));
364 else {
365 G_debug(1, "Cross validation finished correctly");
366
367 G_done_msg(_("Cross validation finished for ew_step = %f and ns_step = %f"), stepE, stepN);
368
369 Segment_close(&in_seg); /* release memory */
370
371 exit(EXIT_SUCCESS);
372 }
373 }
374
375 /* Alloc and load masking matrix */
376 /* encoding: 0 = do not interpolate, 1 = interpolate */
377 have_mask = mask_opt->answer || null_flag->answer;
378 if (have_mask) {
379 int maskfd;
380 int null_count = 0;
381 CELL cval;
382 CELL *maskbuf;
383 char mask_val;
384
385 G_message(_("Mark cells for interpolation"));
386
387 /* use destination window */
388
389 mask_file = G_tempfile();
390 if (Segment_open(&mask_seg, mask_file, nrows, ncols, SEGSIZE, SEGSIZE,
391 sizeof(char), segments_in_memory) != 1)
392 G_fatal_error(_("Can not create temporary file"));
393
394 if (mask_opt->answer) {
395 maskfd = Rast_open_old(mask_opt->answer, "");
396 maskbuf = Rast_allocate_buf(CELL_TYPE);
397 }
398 else {
399 maskfd = -1;
400 maskbuf = NULL;
401 }
402
403 if (null_flag->answer) {
404 inrastfd = Rast_open_old(in_opt->answer, "");
405 drastbuf = Rast_allocate_buf(DCELL_TYPE);
406 }
407 else {
408 inrastfd = -1;
409 drastbuf = NULL;
410 }
411
412 for (row = 0; row < nrows; row++) {
413 G_percent(row, nrows, 9);
414
415 if (mask_opt->answer)
416 Rast_get_c_row(maskfd, maskbuf, row);
417
418 if (null_flag->answer)
419 Rast_get_d_row(inrastfd, drastbuf, row);
420
421 for (col = 0; col < ncols; col++) {
422 mask_val = 1;
423
424 if (mask_opt->answer) {
425 cval = maskbuf[col];
426 if (Rast_is_c_null_value(&cval) || cval == 0)
427 mask_val = 0;
428 }
429
430 if (null_flag->answer && mask_val == 1) {
431 dval = drastbuf[col];
432 if (!Rast_is_d_null_value(&dval))
433 mask_val = 0;
434 else
435 null_count++;
436 }
437 if (Segment_put(&mask_seg, &mask_val, row, col) < 0)
438 G_fatal_error(_("Can not write to temp file"));;
439 }
440 }
441
442 G_percent(row, nrows, 2);
443 if (null_flag->answer) {
444 G_free(drastbuf);
445 Rast_close(inrastfd);
446 }
447 if (mask_opt->answer) {
448 G_free(maskbuf);
449 Rast_close(maskfd);
450 }
451 if (null_flag->answer && null_count == 0 && !mask_opt->answer) {
452 G_fatal_error(_("No NULL cells found in input raster."));
453 }
454 }
455
456 out_file = G_tempfile();
457 if (Segment_open(&out_seg, out_file, nrows, ncols, SEGSIZE, SEGSIZE,
458 sizeof(double), segments_in_memory) != 1)
459 G_fatal_error(_("Can not create temporary file"));
460
461 /* initialize output */
462 G_message(_("Initializing output..."));
463
464 Rast_set_d_null_value(&dval, 1);
465 for (row = 0; row < nrows; row++) {
466 G_percent(row, nrows, 2);
467 for (col = 0; col < ncols; col++) {
468 if (Segment_put(&out_seg, &dval, row, col) < 0)
469 G_fatal_error(_("Can not write to temp file"));;
470 }
471 }
472 G_percent(row, nrows, 2);
473
474 if (grid_opt->answer) {
475 if (0 > Vect_open_new(&Grid, grid_opt->answer, WITH_Z))
476 G_fatal_error(_("Unable to create vector map <%s>"),
477 grid_opt->answer);
478
479 Points = Vect_new_line_struct();
480 Cats = Vect_new_cats_struct();
481 }
482
483 subregion_row = 0;
484 elaboration_reg.south = dest_reg.north;
485 last_row = FALSE;
486 overlap_box.S = dest_box.N;
487 general_box.S = dest_box.N;
488
489 while (last_row == FALSE) { /* For each subregion row */
490 subregion_row++;
491 last_overlap_box.S = overlap_box.S;
492 last_general_box.S = general_box.S;
493 P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
494 GENERAL_ROW);
495
496 /* only works if source reg = dest reg with buffer */
497 /* messing with elaboration region is dangerous... */
498 /* align_elaboration_box(&elaboration_reg, &src_reg, GENERAL_ROW); */
499 align_interp_boxes(&general_box, &overlap_box, &dest_reg,
500 last_general_box, last_overlap_box, GENERAL_ROW);
501
502 if (elaboration_reg.north > dest_reg.north) { /* First row */
503
504 P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
505 FIRST_ROW);
506 /* align_elaboration_box(&elaboration_reg, &src_reg, GENERAL_ROW); */
507 align_interp_boxes(&general_box, &overlap_box, &dest_reg,
508 last_general_box, last_overlap_box, FIRST_ROW);
509 }
510
511 if (elaboration_reg.south <= dest_reg.south) { /* Last row */
512
513 P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
514 LAST_ROW);
515 last_row = TRUE;
516 }
517
518 nsply =
519 ceil((elaboration_reg.north -
520 elaboration_reg.south) / stepN) + 0.5;
521 G_debug(1, "Interpolation: nsply = %d", nsply);
522
523 elaboration_reg.east = dest_reg.west;
524 last_column = FALSE;
525 subregion_col = 0;
526
527 overlap_box.E = dest_box.W;
528 general_box.E = dest_box.W;
529
530 while (last_column == FALSE) { /* For each subregion column */
531 int npoints = 0;
532 int npoints_marked = 1; /* default: all points in output */
533
534 subregion_col++;
535 subregion++;
536 if (nsubregions > 1)
537 G_message(_("subregion %d of %d"), subregion, nsubregions);
538
539 last_overlap_box.E = overlap_box.E;
540 last_general_box.E = general_box.E;
541
542 P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
543 GENERAL_COLUMN);
544
545 /* only works if source reg = dest reg with buffer */
546 /* messing with elaboration region is dangerous... */
547 /* align_elaboration_box(&elaboration_reg, &src_reg, GENERAL_COLUMN); */
548 align_interp_boxes(&general_box, &overlap_box, &dest_reg,
549 last_general_box, last_overlap_box, GENERAL_COLUMN);
550
551 if (elaboration_reg.west < dest_reg.west) { /* First column */
552
553 P_set_regions(&elaboration_reg, &general_box, &overlap_box,
554 dims, FIRST_COLUMN);
555 /* align_elaboration_box(&elaboration_reg, &src_reg, GENERAL_COLUMN); */
556 align_interp_boxes(&general_box, &overlap_box, &dest_reg,
557 last_general_box, last_overlap_box, FIRST_COLUMN);
558 }
559
560 if (elaboration_reg.east >= dest_reg.east) { /* Last column */
561
562 P_set_regions(&elaboration_reg, &general_box, &overlap_box,
563 dims, LAST_COLUMN);
564 last_column = TRUE;
565 }
566 nsplx =
567 ceil((elaboration_reg.east -
568 elaboration_reg.west) / stepE) + 0.5;
569 G_debug(1, "Interpolation: nsplx = %d", nsplx);
570
571 if (grid_opt->answer) {
572 Vect_reset_cats(Cats);
573 Vect_cat_set(Cats, 1, cat++);
574 Vect_reset_line(Points);
575 Vect_append_point(Points, general_box.W, general_box.S, 0);
576 Vect_append_point(Points, general_box.E, general_box.S, 0);
577 Vect_append_point(Points, general_box.E, general_box.N, 0);
578 Vect_append_point(Points, general_box.W, general_box.N, 0);
579 Vect_append_point(Points, general_box.W, general_box.S, 0);
580 Vect_write_line(&Grid, GV_LINE, Points, Cats);
581 Vect_reset_line(Points);
582 Vect_append_point(Points, (general_box.E + general_box.W) / 2,
583 (general_box.N + general_box.S) / 2, 0);
584 Vect_write_line(&Grid, GV_POINT, Points, Cats);
585 }
586
587 /* reading points in interpolation region */
588 G_debug(1, "reading points from input raster...");
589 dim_vect = nsplx * nsply;
590
591 observ =
592 P_Read_Raster_Region_Map(&in_seg, &elaboration_reg,
593 &src_reg, &npoints, dim_vect);
594
595 G_debug(1, "%d valid points", npoints);
596
597 G_debug(1,
598 "Interpolation: (%d,%d): Number of points in <elaboration_box> is %d",
599 subregion_row, subregion_col, npoints);
600
601 /* Mean calculation for observed non-NULL points */
602 if (npoints > 0)
603 mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
604 else
605 mean = 0;
606 G_debug(1, "Interpolation: (%d,%d): mean=%lf",
607 subregion_row, subregion_col, mean);
608
609 observ_marked = NULL;
610
611 if (have_mask) {
612 /* collect unmasked output cells */
613
614 G_debug(1, "collect unmasked output cells");
615
616 observ_marked =
617 P_Read_Raster_Region_masked(&mask_seg, &dest_reg,
618 dest_box, general_box,
619 &npoints_marked, dim_vect, mean);
620
621 G_debug(1, "%d cells marked in general", npoints_marked);
622 if (npoints_marked == 0) {
623 G_free(observ_marked);
624 observ_marked = NULL;
625 npoints = 1; /* disable warning below */
626 }
627 }
628
629 if (npoints > 0 && npoints_marked > 0) { /* */
630 int i;
631
632 nparameters = nsplx * nsply;
633 BW = P_get_BandWidth(interp_method, nsply > nsplx ? nsply : nsplx);
634
635 /* Least Squares system */
636 N = G_alloc_matrix(nparameters, BW); /* Normal matrix */
637 TN = G_alloc_vector(nparameters); /* vector */
638 parVect = G_alloc_vector(nparameters); /* Parameters vector */
639 obsVect = G_alloc_matrix(npoints, 3); /* Observation vector */
640 Q = G_alloc_vector(npoints); /* "a priori" var-cov matrix */
641
642 for (i = 0; i < npoints; i++) { /* Setting obsVect vector & Q matrix */
643 Q[i] = 1; /* Q=I */
644 obsVect[i][0] = observ[i].coordX;
645 obsVect[i][1] = observ[i].coordY;
646 obsVect[i][2] = observ[i].coordZ - mean;
647 }
648
649 G_free(observ);
650
651 /* Bilinear interpolation */
652 if (interp_method == P_BILINEAR) {
653 G_debug(1,
654 "Interpolation: (%d,%d): Bilinear interpolation...",
655 subregion_row, subregion_col);
656 normalDefBilin(N, TN, Q, obsVect, stepE, stepN, nsplx,
657 nsply, elaboration_reg.west,
658 elaboration_reg.south, npoints,
659 nparameters, BW);
660 nCorrectGrad(N, lambda, nsplx, nsply, stepE, stepN);
661 }
662 /* Bicubic interpolation */
663 else {
664 G_debug(1,
665 "Interpolation: (%d,%d): Bicubic interpolation...",
666 subregion_row, subregion_col);
667 normalDefBicubic(N, TN, Q, obsVect, stepE, stepN, nsplx,
668 nsply, elaboration_reg.west,
669 elaboration_reg.south, npoints,
670 nparameters, BW);
671 nCorrectGrad(N, lambda, nsplx, nsply, stepE, stepN);
672 }
673
674 G_math_solver_cholesky_sband(N, parVect, TN, nparameters, BW);
675
676 G_free_matrix(N);
677 G_free_vector(TN);
678 G_free_vector(Q);
679
680 if (!observ_marked) { /* interpolate full output raster */
681 G_debug(1, "Interpolation: (%d,%d): Regular_Points...",
682 subregion_row, subregion_col);
683
684 P_Regular_Points(&elaboration_reg, &dest_reg, general_box,
685 overlap_box, &out_seg,
686 parVect, stepN, stepE, dims.overlap, mean,
687 nsplx, nsply, nrows, ncols, interp_method);
688 }
689 else { /* only interpolate selected cells */
690
691 G_debug(1, "Interpolation of %d selected cells...",
692 npoints_marked);
693
694 P_Sparse_Raster_Points(&out_seg,
695 &elaboration_reg, &dest_reg,
696 general_box, overlap_box,
697 observ_marked, parVect,
698 stepE, stepN,
699 dims.overlap, nsplx, nsply,
700 npoints_marked, interp_method, mean);
701
702 G_free(observ_marked);
703 } /* end NULL cells */
704 G_free_vector(parVect);
705 G_free_matrix(obsVect);
706 }
707 else {
708 if (observ)
709 G_free(observ);
710 if (npoints == 0)
711 G_warning(_("No data within this subregion. "
712 "Consider increasing the spline step."));
713 }
714 } /*! END WHILE; last_column = TRUE */
715 } /*! END WHILE; last_row = TRUE */
716
717 Segment_close(&in_seg); /* release memory */
718
719 if (have_mask) {
720 Segment_close(&mask_seg); /* release memory */
721 }
722
723 G_message(_("Writing output..."));
724 /* Writing the output raster map */
725 Rast_set_fp_type(DCELL_TYPE);
726 outrastfd = Rast_open_fp_new(out_opt->answer);
727
728 /* check */
729 {
730 int nonulls = 0;
731
732 drastbuf = Rast_allocate_d_buf();
733
734 for (row = 0; row < dest_reg.rows; row++) {
735 G_percent(row, dest_reg.rows, 9);
736 for (col = 0; col < dest_reg.cols; col++) {
737 Segment_get(&out_seg, &dval, row, col);
738 drastbuf[col] = dval;
739 if (!Rast_is_d_null_value(&dval))
740 nonulls++;
741 }
742 Rast_put_d_row(outrastfd, drastbuf);
743 }
744 G_percent(1, 1, 1);
745 if (!nonulls)
746 G_warning("only NULL cells in output raster");
747 }
748
749 Segment_close(&out_seg); /* release memory */
750
751 Rast_close(outrastfd);
752
753 /* set map title */
754 sprintf(title, "%s interpolation with Tykhonov regularization",
755 method_opt->answer);
756 Rast_put_cell_title(out_opt->answer, title);
757 /* write map history */
758 Rast_short_history(out_opt->answer, "raster", &history);
759 Rast_command_history(&history);
760 Rast_write_history(out_opt->answer, &history);
761
762 if (grid_opt->answer) {
763 Vect_build(&Grid);
764 Vect_hist_command(&Grid);
765 Vect_close(&Grid);
766 }
767
768 G_done_msg(" ");
769
770 exit(EXIT_SUCCESS);
771 } /*END MAIN */
772