source: grass/trunk/vector/v.surf.bspline/main.c

Last change on this file was 73507, checked in by mmetz, 6 years ago

v.surf.bspline: use Segment_[open|close]()
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id
  • Property svn:mime-type set to text/x-csrc
File size: 29.1 KB
Line 
1
2/**********************************************************************
3 *
4 * MODULE: v.surf.bspline
5 *
6 * AUTHOR(S): Roberto Antolin & Gonzalo Moreno
7 * update for grass7 by Markus Metz
8 *
9 * PURPOSE: Spline Interpolation
10 *
11 * COPYRIGHT: (C) 2006 by Politecnico di Milano -
12 * Polo Regionale di Como
13 *
14 * This program is free software under the
15 * GNU General Public License (>=v2).
16 * Read the file COPYING that comes with GRASS
17 * for details.
18 *
19 **********************************************************************/
20
21/* INCLUDES */
22#include <stdlib.h>
23#include <string.h>
24#include <math.h>
25#include <sys/types.h>
26#include <sys/stat.h>
27#include <fcntl.h>
28#include "bspline.h"
29#include <grass/N_pde.h>
30
31#define SEGSIZE 64
32
33/* GLOBAL VARIABLES */
34int bspline_field;
35char *bspline_column;
36
37/*--------------------------------------------------------------------*/
38int main(int argc, char *argv[])
39{
40 /* Variable declarations */
41 int nsply, nsplx, nrows, ncols, nsplx_adj, nsply_adj;
42 int nsubregion_col, nsubregion_row, subregion_row, subregion_col;
43 int subregion = 0, nsubregions = 0;
44 int last_row, last_column, grid, bilin, ext, flag_auxiliar, cross; /* booleans */
45 double stepN, stepE, lambda, mean;
46 double N_extension, E_extension, edgeE, edgeN;
47
48 const char *mapset, *drv, *db, *vector, *map;
49 char table_name[GNAME_MAX], title[64];
50 char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
51
52 int dim_vect, nparameters, BW;
53 int *lineVect; /* Vector restoring primitive's ID */
54 double *TN, *Q, *parVect; /* Interpolating and least-square vectors */
55 double **N, **obsVect; /* Interpolation and least-square matrix */
56
57 SEGMENT out_seg, mask_seg;
58 char *out_file, *mask_file;
59 double seg_size;
60 int seg_mb, segments_in_memory;
61 int have_mask;
62
63 /* Structs declarations */
64 int raster;
65 struct Map_info In, In_ext, Out;
66 struct History history;
67
68 struct GModule *module;
69 struct Option *in_opt, *in_ext_opt, *out_opt, *out_map_opt, *stepE_opt,
70 *stepN_opt, *lambda_f_opt, *type_opt, *dfield_opt, *col_opt, *mask_opt,
71 *memory_opt, *solver, *error, *iter;
72 struct Flag *cross_corr_flag, *spline_step_flag;
73
74 struct Reg_dimens dims;
75 struct Cell_head elaboration_reg, original_reg;
76 struct bound_box general_box, overlap_box, original_box;
77
78 struct Point *observ;
79 struct line_cats *Cats;
80 dbCatValArray cvarr;
81
82 int with_z;
83 int nrec, ctype = 0;
84 struct field_info *Fi;
85 dbDriver *driver, *driver_cats;
86
87 /*----------------------------------------------------------------*/
88 /* Options declarations */
89 module = G_define_module();
90 G_add_keyword(_("vector"));
91 G_add_keyword(_("surface"));
92 G_add_keyword(_("interpolation"));
93 G_add_keyword(_("LIDAR"));
94 module->description =
95 _("Performs bicubic or bilinear spline interpolation with Tykhonov regularization.");
96
97 cross_corr_flag = G_define_flag();
98 cross_corr_flag->key = 'c';
99 cross_corr_flag->description =
100 _("Find the best Tykhonov regularizing parameter using a \"leave-one-out\" cross validation method");
101
102 spline_step_flag = G_define_flag();
103 spline_step_flag->key = 'e';
104 spline_step_flag->label = _("Estimate point density and distance");
105 spline_step_flag->description =
106 _("Estimate point density and distance in map units for the input vector points within the current region extents and quit");
107
108 in_opt = G_define_standard_option(G_OPT_V_INPUT);
109 in_opt->label = _("Name of input vector point map");
110
111 dfield_opt = G_define_standard_option(G_OPT_V_FIELD);
112 dfield_opt->guisection = _("Settings");
113
114 col_opt = G_define_standard_option(G_OPT_DB_COLUMN);
115 col_opt->required = NO;
116 col_opt->label =
117 _("Name of the attribute column with values to be used for approximation");
118 col_opt->description = _("If not given and input is 3D vector map then z-coordinates are used.");
119 col_opt->guisection = _("Settings");
120
121 in_ext_opt = G_define_standard_option(G_OPT_V_INPUT);
122 in_ext_opt->key = "sparse_input";
123 in_ext_opt->required = NO;
124 in_ext_opt->label =
125 _("Name of input vector map with sparse points");
126
127 out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
128 out_opt->required = NO;
129 out_opt->guisection = _("Outputs");
130
131 out_map_opt = G_define_standard_option(G_OPT_R_OUTPUT);
132 out_map_opt->key = "raster_output";
133 out_map_opt->required = NO;
134 out_map_opt->guisection = _("Outputs");
135
136 mask_opt = G_define_standard_option(G_OPT_R_INPUT);
137 mask_opt->key = "mask";
138 mask_opt->label = _("Raster map to use for masking (applies to raster output only)");
139 mask_opt->description = _("Only cells that are not NULL and not zero are interpolated");
140 mask_opt->required = NO;
141
142 stepE_opt = G_define_option();
143 stepE_opt->key = "ew_step";
144 stepE_opt->type = TYPE_DOUBLE;
145 stepE_opt->required = NO;
146 stepE_opt->label =
147 _("Length of each spline step in the east-west direction");
148 stepE_opt->description = _("Default: 4 * east-west resolution");
149 stepE_opt->guisection = _("Settings");
150
151 stepN_opt = G_define_option();
152 stepN_opt->key = "ns_step";
153 stepN_opt->type = TYPE_DOUBLE;
154 stepN_opt->required = NO;
155 stepN_opt->label =
156 _("Length of each spline step in the north-south direction");
157 stepN_opt->description = _("Default: 4 * north-south resolution");
158 stepN_opt->guisection = _("Settings");
159
160 type_opt = G_define_option();
161 type_opt->key = "method";
162 type_opt->description = _("Spline interpolation algorithm");
163 type_opt->type = TYPE_STRING;
164 type_opt->options = "bilinear,bicubic";
165 type_opt->answer = "bilinear";
166 type_opt->guisection = _("Settings");
167 G_asprintf((char **) &(type_opt->descriptions),
168 "bilinear;%s;bicubic;%s",
169 _("Bilinear interpolation"),
170 _("Bicubic interpolation"));
171
172 lambda_f_opt = G_define_option();
173 lambda_f_opt->key = "lambda_i";
174 lambda_f_opt->type = TYPE_DOUBLE;
175 lambda_f_opt->required = NO;
176 lambda_f_opt->description = _("Tykhonov regularization parameter (affects smoothing)");
177 lambda_f_opt->answer = "0.01";
178 lambda_f_opt->guisection = _("Settings");
179
180 solver = N_define_standard_option(N_OPT_SOLVER_SYMM);
181 solver->options = "cholesky,cg";
182 solver->answer = "cholesky";
183
184 iter = N_define_standard_option(N_OPT_MAX_ITERATIONS);
185
186 error = N_define_standard_option(N_OPT_ITERATION_ERROR);
187
188 memory_opt = G_define_option();
189 memory_opt->key = "memory";
190 memory_opt->type = TYPE_INTEGER;
191 memory_opt->required = NO;
192 memory_opt->answer = "300";
193 memory_opt->label = _("Maximum memory to be used (in MB)");
194 memory_opt->description = _("Cache size for raster rows");
195
196 /*----------------------------------------------------------------*/
197 /* Parsing */
198 G_gisinit(argv[0]);
199 if (G_parser(argc, argv))
200 exit(EXIT_FAILURE);
201
202 vector = out_opt->answer;
203 map = out_map_opt->answer;
204
205 if (vector && map)
206 G_fatal_error(_("Choose either vector or raster output, not both"));
207
208 if (!vector && !map && !cross_corr_flag->answer)
209 G_fatal_error(_("No raster or vector or cross-validation output"));
210
211 if (!strcmp(type_opt->answer, "linear"))
212 bilin = P_BILINEAR;
213 else
214 bilin = P_BICUBIC;
215
216 G_get_set_window(&original_reg);
217 stepN = 4 * original_reg.ns_res;
218 if (stepN_opt->answer)
219 stepN = atof(stepN_opt->answer);
220 stepE = 4 * original_reg.ew_res;
221 if (stepE_opt->answer)
222 stepE = atof(stepE_opt->answer);
223 lambda = atof(lambda_f_opt->answer);
224
225 flag_auxiliar = FALSE;
226
227 drv = db_get_default_driver_name();
228 if (!drv) {
229 if (db_set_default_connection() != DB_OK)
230 G_fatal_error(_("Unable to set default DB connection"));
231 drv = db_get_default_driver_name();
232 }
233 db = db_get_default_database_name();
234 if (!db)
235 G_fatal_error(_("No default DB defined"));
236
237 /* Set auxiliary table's name */
238 if (vector) {
239 if (G_name_is_fully_qualified(out_opt->answer, xname, xmapset)) {
240 sprintf(table_name, "%s_aux", xname);
241 }
242 else
243 sprintf(table_name, "%s_aux", out_opt->answer);
244 }
245
246 /* Something went wrong in a previous v.surf.bspline execution */
247 if (db_table_exists(drv, db, table_name)) {
248 /* Start driver and open db */
249 driver = db_start_driver_open_database(drv, db);
250 if (driver == NULL)
251 G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
252 drv);
253 db_set_error_handler_driver(driver);
254
255 if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
256 G_fatal_error(_("Old auxiliary table could not be dropped"));
257 db_close_database_shutdown_driver(driver);
258 }
259
260 /* Open input vector */
261 if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL)
262 G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
263
264 Vect_set_open_level(1); /* WITHOUT TOPOLOGY */
265 if (1 > Vect_open_old(&In, in_opt->answer, mapset))
266 G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
267 in_opt->answer);
268
269 bspline_field = 0; /* assume 3D input */
270 bspline_column = col_opt->answer;
271
272 with_z = !bspline_column && Vect_is_3d(&In);
273
274 if (Vect_is_3d(&In)) {
275 if (!with_z)
276 G_verbose_message(_("Input is 3D: using attribute values instead of z-coordinates for approximation"));
277 else
278 G_verbose_message(_("Input is 3D: using z-coordinates for approximation"));
279 }
280 else { /* 2D */
281 if (!bspline_column)
282 G_fatal_error(_("Input vector map is 2D. Parameter <%s> required."), col_opt->key);
283 }
284
285 if (!with_z) {
286 bspline_field = Vect_get_field_number(&In, dfield_opt->answer);
287 }
288
289 /* Estimate point density and mean distance for current region */
290 if (spline_step_flag->answer) {
291 double dens, dist;
292
293 if (P_estimate_splinestep(&In, &dens, &dist) == 0) {
294 fprintf(stdout, _("Estimated point density: %.4g"), dens);
295 fprintf(stdout, _("Estimated mean distance between points: %.4g"), dist);
296 }
297 else {
298 fprintf(stdout, _("No points in current region"));
299 }
300
301 Vect_close(&In);
302 exit(EXIT_SUCCESS);
303 }
304
305 /*----------------------------------------------------------------*/
306 /* Cross-correlation begins */
307 if (cross_corr_flag->answer) {
308 G_debug(1, "CrossCorrelation()");
309 cross = cross_correlation(&In, stepE, stepN);
310
311 if (cross != TRUE)
312 G_fatal_error(_("Cross validation didn't finish correctly"));
313 else {
314 G_debug(1, "Cross validation finished correctly");
315
316 Vect_close(&In);
317
318 G_done_msg(_("Cross validation finished for ew_step = %f and ns_step = %f"), stepE, stepN);
319 exit(EXIT_SUCCESS);
320 }
321 }
322
323 /* Open input ext vector */
324 ext = FALSE;
325 if (in_ext_opt->answer) {
326 ext = TRUE;
327 G_message(_("Vector map <%s> of sparse points will be interpolated"),
328 in_ext_opt->answer);
329
330 if ((mapset = G_find_vector2(in_ext_opt->answer, "")) == NULL)
331 G_fatal_error(_("Vector map <%s> not found"), in_ext_opt->answer);
332
333 Vect_set_open_level(1); /* WITHOUT TOPOLOGY */
334 if (1 > Vect_open_old(&In_ext, in_ext_opt->answer, mapset))
335 G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
336 in_opt->answer);
337 }
338
339 /* Open output map */
340 /* vector output */
341 if (vector && !map) {
342 if (strcmp(drv, "dbf") == 0)
343 G_fatal_error(_("Sorry, the <%s> driver is not compatible with "
344 "the vector output of this module. "
345 "Try with raster output or another driver."), drv);
346
347 Vect_check_input_output_name(in_opt->answer, out_opt->answer,
348 G_FATAL_EXIT);
349 grid = FALSE;
350
351 if (0 > Vect_open_new(&Out, out_opt->answer, WITH_Z))
352 G_fatal_error(_("Unable to create vector map <%s>"),
353 out_opt->answer);
354
355 /* Copy vector Head File */
356 if (ext == FALSE) {
357 Vect_copy_head_data(&In, &Out);
358 Vect_hist_copy(&In, &Out);
359 }
360 else {
361 Vect_copy_head_data(&In_ext, &Out);
362 Vect_hist_copy(&In_ext, &Out);
363 }
364 Vect_hist_command(&Out);
365
366 G_verbose_message(_("Points in input vector map <%s> will be interpolated"),
367 vector);
368 }
369
370
371 /* read z values from attribute table */
372 if (bspline_field > 0) {
373 G_message(_("Reading values from attribute table..."));
374 db_CatValArray_init(&cvarr);
375 Fi = Vect_get_field(&In, bspline_field);
376 if (Fi == NULL)
377 G_fatal_error(_("Cannot read layer info"));
378
379 driver_cats = db_start_driver_open_database(Fi->driver, Fi->database);
380 /*G_debug (0, _("driver=%s db=%s"), Fi->driver, Fi->database); */
381
382 if (driver_cats == NULL)
383 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
384 Fi->database, Fi->driver);
385 db_set_error_handler_driver(driver_cats);
386
387 nrec =
388 db_select_CatValArray(driver_cats, Fi->table, Fi->key,
389 col_opt->answer, NULL, &cvarr);
390 G_debug(3, "nrec = %d", nrec);
391
392 ctype = cvarr.ctype;
393 if (ctype != DB_C_TYPE_INT && ctype != DB_C_TYPE_DOUBLE)
394 G_fatal_error(_("Column type not supported"));
395
396 if (nrec < 0)
397 G_fatal_error(_("Unable to select data from table"));
398
399 G_verbose_message(_("%d records selected from table"), nrec);
400
401 db_close_database_shutdown_driver(driver_cats);
402 }
403
404 /*----------------------------------------------------------------*/
405 /* Interpolation begins */
406 G_debug(1, "Interpolation()");
407
408 /* Open driver and database */
409 driver = db_start_driver_open_database(drv, db);
410 if (driver == NULL)
411 G_fatal_error(_("No database connection for driver <%s> is defined. "
412 "Run db.connect."), drv);
413 db_set_error_handler_driver(driver);
414
415 /* Create auxiliary table */
416 if (vector) {
417 if ((flag_auxiliar = P_Create_Aux4_Table(driver, table_name)) == FALSE) {
418 P_Drop_Aux_Table(driver, table_name);
419 G_fatal_error(_("Interpolation: Creating table: "
420 "It was impossible to create table <%s>."),
421 table_name);
422 }
423 /* db_create_index2(driver, table_name, "ID"); */
424 /* sqlite likes that ??? */
425 db_close_database_shutdown_driver(driver);
426 driver = db_start_driver_open_database(drv, db);
427 }
428
429 /* raster output */
430 raster = -1;
431 Rast_set_fp_type(DCELL_TYPE);
432 if (!vector && map) {
433 grid = TRUE;
434 raster = Rast_open_fp_new(out_map_opt->answer);
435
436 G_verbose_message(_("Cells for raster map <%s> will be interpolated"),
437 map);
438 }
439
440 /* Setting regions and boxes */
441 G_debug(1, "Interpolation: Setting regions and boxes");
442 G_get_window(&elaboration_reg);
443 Vect_region_box(&original_reg, &original_box);
444 Vect_region_box(&elaboration_reg, &overlap_box);
445 Vect_region_box(&elaboration_reg, &general_box);
446
447 nrows = Rast_window_rows();
448 ncols = Rast_window_cols();
449
450 /* Alloc raster matrix */
451 have_mask = 0;
452 out_file = mask_file = NULL;
453 if (grid == TRUE) {
454 int row;
455 DCELL *drastbuf;
456
457 seg_mb = atoi(memory_opt->answer);
458 if (seg_mb < 3)
459 G_fatal_error(_("Memory in MB must be >= 3"));
460
461 if (mask_opt->answer)
462 seg_size = sizeof(double) + sizeof(char);
463 else
464 seg_size = sizeof(double);
465
466 seg_size = (seg_size * SEGSIZE * SEGSIZE) / (1 << 20);
467 segments_in_memory = seg_mb / seg_size + 0.5;
468 G_debug(1, "%d %dx%d segments held in memory", segments_in_memory, SEGSIZE, SEGSIZE);
469
470 out_file = G_tempfile();
471 if (Segment_open(&out_seg, out_file, nrows, ncols, SEGSIZE,
472 SEGSIZE, sizeof(double), segments_in_memory) != 1)
473 G_fatal_error(_("Can not create temporary file"));
474
475 /* initialize output */
476 G_message(_("Initializing output..."));
477
478 drastbuf = Rast_allocate_buf(DCELL_TYPE);
479 Rast_set_d_null_value(drastbuf, ncols);
480 for (row = 0; row < nrows; row++) {
481 G_percent(row, nrows, 2);
482 Segment_put_row(&out_seg, drastbuf, row);
483 }
484 G_percent(row, nrows, 2);
485
486 if (mask_opt->answer) {
487 int col, maskfd;
488 DCELL dval;
489 char mask_val;
490
491 G_message(_("Load masking map"));
492
493 mask_file = G_tempfile();
494 if (Segment_open(&mask_seg, mask_file, nrows, ncols, SEGSIZE,
495 SEGSIZE, sizeof(char), segments_in_memory) != 1)
496 G_fatal_error(_("Can not create temporary file"));
497
498 maskfd = Rast_open_old(mask_opt->answer, "");
499
500 for (row = 0; row < nrows; row++) {
501 G_percent(row, nrows, 2);
502 Rast_get_d_row(maskfd, drastbuf, row);
503 for (col = 0; col < ncols; col++) {
504 dval = drastbuf[col];
505 if (Rast_is_d_null_value(&dval) || dval == 0)
506 mask_val = 0;
507 else
508 mask_val = 1;
509
510 Segment_put(&mask_seg, &mask_val, row, col);
511 }
512 }
513
514 G_percent(row, nrows, 2);
515 Rast_close(maskfd);
516
517 have_mask = 1;
518 }
519 G_free(drastbuf);
520 }
521
522 /*------------------------------------------------------------------
523 | Subdividing and working with tiles:
524 | Each original region will be divided into several subregions.
525 | Each one will be overlapped by its neighbouring subregions.
526 | The overlapping is calculated as a fixed OVERLAP_SIZE times
527 | the largest spline step plus 2 * edge
528 ----------------------------------------------------------------*/
529
530 /* Fixing parameters of the elaboration region */
531 P_zero_dim(&dims); /* Set dim struct to zero */
532
533 nsplx_adj = NSPLX_MAX;
534 nsply_adj = NSPLY_MAX;
535 if (stepN > stepE)
536 dims.overlap = OVERLAP_SIZE * stepN;
537 else
538 dims.overlap = OVERLAP_SIZE * stepE;
539 P_get_edge(bilin, &dims, stepE, stepN);
540 P_set_dim(&dims, stepE, stepN, &nsplx_adj, &nsply_adj);
541
542 G_verbose_message(_("Adjusted EW splines %d"), nsplx_adj);
543 G_verbose_message(_("Adjusted NS splines %d"), nsply_adj);
544
545 /* calculate number of subregions */
546 edgeE = dims.ew_size - dims.overlap - 2 * dims.edge_v;
547 edgeN = dims.sn_size - dims.overlap - 2 * dims.edge_h;
548
549 N_extension = original_reg.north - original_reg.south;
550 E_extension = original_reg.east - original_reg.west;
551
552 nsubregion_col = ceil(E_extension / edgeE) + 0.5;
553 nsubregion_row = ceil(N_extension / edgeN) + 0.5;
554
555 if (nsubregion_col < 0)
556 nsubregion_col = 0;
557 if (nsubregion_row < 0)
558 nsubregion_row = 0;
559
560 nsubregions = nsubregion_row * nsubregion_col;
561
562 /* Creating line and categories structs */
563 Cats = Vect_new_cats_struct();
564 Vect_cat_set(Cats, 1, 0);
565
566 subregion_row = 0;
567 elaboration_reg.south = original_reg.north;
568 last_row = FALSE;
569
570 while (last_row == FALSE) { /* For each subregion row */
571 subregion_row++;
572 P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
573 GENERAL_ROW);
574
575 if (elaboration_reg.north > original_reg.north) { /* First row */
576
577 P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
578 FIRST_ROW);
579 }
580
581 if (elaboration_reg.south <= original_reg.south) { /* Last row */
582
583 P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
584 LAST_ROW);
585 last_row = TRUE;
586 }
587
588 nsply =
589 ceil((elaboration_reg.north -
590 elaboration_reg.south) / stepN) + 0.5;
591 G_debug(1, "Interpolation: nsply = %d", nsply);
592 /*
593 if (nsply > NSPLY_MAX)
594 nsply = NSPLY_MAX;
595 */
596 elaboration_reg.east = original_reg.west;
597 last_column = FALSE;
598 subregion_col = 0;
599
600 /* TODO: process each subregion using its own thread (via OpenMP or pthreads) */
601 /* I'm not sure about pthreads, but you can tell OpenMP to start all at the
602 same time and it will keep num_workers supplied with the next job as free
603 cpus become available */
604 while (last_column == FALSE) { /* For each subregion column */
605 int npoints = 0;
606 /* needed for sparse points interpolation */
607 int npoints_ext, *lineVect_ext = NULL;
608 double **obsVect_ext; /*, mean_ext = .0; */
609 struct Point *observ_ext;
610
611 subregion_col++;
612 subregion++;
613 if (nsubregions > 1)
614 G_message(_("Processing subregion %d of %d..."), subregion, nsubregions);
615
616 P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
617 GENERAL_COLUMN);
618
619 if (elaboration_reg.west < original_reg.west) { /* First column */
620
621 P_set_regions(&elaboration_reg, &general_box, &overlap_box,
622 dims, FIRST_COLUMN);
623 }
624
625 if (elaboration_reg.east >= original_reg.east) { /* Last column */
626
627 P_set_regions(&elaboration_reg, &general_box, &overlap_box,
628 dims, LAST_COLUMN);
629 last_column = TRUE;
630 }
631 nsplx =
632 ceil((elaboration_reg.east -
633 elaboration_reg.west) / stepE) + 0.5;
634 G_debug(1, "Interpolation: nsplx = %d", nsplx);
635 /*
636 if (nsplx > NSPLX_MAX)
637 nsplx = NSPLX_MAX;
638 */
639 G_debug(1, "Interpolation: (%d,%d): subregion bounds",
640 subregion_row, subregion_col);
641 G_debug(1, "Interpolation: \t\tNORTH:%.2f\t",
642 elaboration_reg.north);
643 G_debug(1, "Interpolation: WEST:%.2f\t\tEAST:%.2f",
644 elaboration_reg.west, elaboration_reg.east);
645 G_debug(1, "Interpolation: \t\tSOUTH:%.2f",
646 elaboration_reg.south);
647
648#ifdef DEBUG_SUBREGIONS
649 fprintf(stdout, "B 5\n");
650 fprintf(stdout, " %.11g %.11g\n", elaboration_reg.east, elaboration_reg.north);
651 fprintf(stdout, " %.11g %.11g\n", elaboration_reg.west, elaboration_reg.north);
652 fprintf(stdout, " %.11g %.11g\n", elaboration_reg.west, elaboration_reg.south);
653 fprintf(stdout, " %.11g %.11g\n", elaboration_reg.east, elaboration_reg.south);
654 fprintf(stdout, " %.11g %.11g\n", elaboration_reg.east, elaboration_reg.north);
655 fprintf(stdout, "C 1 1\n");
656 fprintf(stdout, " %.11g %.11g\n", (elaboration_reg.west + elaboration_reg.east) / 2,
657 (elaboration_reg.south + elaboration_reg.north) / 2);
658 fprintf(stdout, " 1 %d\n", subregion);
659#endif
660
661
662
663 /* reading points in interpolation region */
664 dim_vect = nsplx * nsply;
665 observ_ext = NULL;
666 if (grid == FALSE && ext == TRUE) {
667 observ_ext =
668 P_Read_Vector_Region_Map(&In_ext,
669 &elaboration_reg,
670 &npoints_ext, dim_vect,
671 1);
672 }
673 else
674 npoints_ext = 1;
675
676 if (grid == TRUE && have_mask) {
677 /* any unmasked cells in general region ? */
678 mean = 0;
679 observ_ext =
680 P_Read_Raster_Region_masked(&mask_seg, &original_reg,
681 original_box, general_box,
682 &npoints_ext, dim_vect, mean);
683 }
684
685 observ = NULL;
686 if (npoints_ext > 0) {
687 observ =
688 P_Read_Vector_Region_Map(&In, &elaboration_reg, &npoints,
689 dim_vect, bspline_field);
690 }
691 else
692 npoints = 1;
693
694 G_debug(1,
695 "Interpolation: (%d,%d): Number of points in <elaboration_box> is %d",
696 subregion_row, subregion_col, npoints);
697 if (npoints > 0)
698 G_verbose_message(_("%d points found in this subregion"), npoints);
699 /* only interpolate if there are any points in current subregion */
700 if (npoints > 0 && npoints_ext > 0) {
701 int i;
702
703 nparameters = nsplx * nsply;
704 BW = P_get_BandWidth(bilin, nsply);
705
706 /* Least Squares system */
707 N = G_alloc_matrix(nparameters, BW); /* Normal matrix */
708 TN = G_alloc_vector(nparameters); /* vector */
709 parVect = G_alloc_vector(nparameters); /* Parameters vector */
710 obsVect = G_alloc_matrix(npoints, 3); /* Observation vector */
711 Q = G_alloc_vector(npoints); /* "a priori" var-cov matrix */
712 lineVect = G_alloc_ivector(npoints); /* */
713
714 for (i = 0; i < npoints; i++) { /* Setting obsVect vector & Q matrix */
715 double dval;
716
717 Q[i] = 1; /* Q=I */
718 lineVect[i] = observ[i].lineID;
719 obsVect[i][0] = observ[i].coordX;
720 obsVect[i][1] = observ[i].coordY;
721
722 /* read z coordinates from attribute table */
723 if (bspline_field > 0) {
724 int cat, ival, ret;
725
726 cat = observ[i].cat;
727 if (cat < 0)
728 continue;
729
730 if (ctype == DB_C_TYPE_INT) {
731 ret =
732 db_CatValArray_get_value_int(&cvarr, cat,
733 &ival);
734 obsVect[i][2] = ival;
735 observ[i].coordZ = ival;
736 }
737 else { /* DB_C_TYPE_DOUBLE */
738 ret =
739 db_CatValArray_get_value_double(&cvarr, cat,
740 &dval);
741 obsVect[i][2] = dval;
742 observ[i].coordZ = dval;
743 }
744 if (ret != DB_OK) {
745 G_warning(_("Interpolation: (%d,%d): No record for point (cat = %d)"),
746 subregion_row, subregion_col, cat);
747 continue;
748 }
749 }
750 /* use z coordinates of 3D vector */
751 else {
752 obsVect[i][2] = observ[i].coordZ;
753 }
754 }
755
756 /* Mean calculation for every point */
757 mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
758
759 G_debug(1, "Interpolation: (%d,%d): mean=%lf",
760 subregion_row, subregion_col, mean);
761
762 G_free(observ);
763
764 for (i = 0; i < npoints; i++)
765 obsVect[i][2] -= mean;
766
767 /* Bilinear interpolation */
768 if (bilin) {
769 G_debug(1,
770 "Interpolation: (%d,%d): Bilinear interpolation...",
771 subregion_row, subregion_col);
772 normalDefBilin(N, TN, Q, obsVect, stepE, stepN, nsplx,
773 nsply, elaboration_reg.west,
774 elaboration_reg.south, npoints,
775 nparameters, BW);
776 nCorrectGrad(N, lambda, nsplx, nsply, stepE, stepN);
777 }
778 /* Bicubic interpolation */
779 else {
780 G_debug(1,
781 "Interpolation: (%d,%d): Bicubic interpolation...",
782 subregion_row, subregion_col);
783 normalDefBicubic(N, TN, Q, obsVect, stepE, stepN, nsplx,
784 nsply, elaboration_reg.west,
785 elaboration_reg.south, npoints,
786 nparameters, BW);
787 nCorrectGrad(N, lambda, nsplx, nsply, stepE, stepN);
788 }
789
790 if(G_strncasecmp(solver->answer, "cg", 2) == 0)
791 G_math_solver_cg_sband(N, parVect, TN, nparameters, BW, atoi(iter->answer), atof(error->answer));
792 else
793 G_math_solver_cholesky_sband(N, parVect, TN, nparameters, BW);
794
795
796 G_free_matrix(N);
797 G_free_vector(TN);
798 G_free_vector(Q);
799
800 if (grid == TRUE) { /* GRID INTERPOLATION ==> INTERPOLATION INTO A RASTER */
801 G_debug(1, "Interpolation: (%d,%d): Regular_Points...",
802 subregion_row, subregion_col);
803
804 if (!have_mask) {
805 P_Regular_Points(&elaboration_reg, &original_reg, general_box,
806 overlap_box, &out_seg, parVect,
807 stepN, stepE, dims.overlap, mean,
808 nsplx, nsply, nrows, ncols, bilin);
809 }
810 else {
811 P_Sparse_Raster_Points(&out_seg,
812 &elaboration_reg, &original_reg,
813 general_box, overlap_box,
814 observ_ext, parVect,
815 stepE, stepN,
816 dims.overlap, nsplx, nsply,
817 npoints_ext, bilin, mean);
818 }
819 }
820 else { /* OBSERVATION POINTS INTERPOLATION */
821 if (ext == FALSE) {
822 G_debug(1, "Interpolation: (%d,%d): Sparse_Points...",
823 subregion_row, subregion_col);
824 P_Sparse_Points(&Out, &elaboration_reg, general_box,
825 overlap_box, obsVect, parVect,
826 lineVect, stepE, stepN,
827 dims.overlap, nsplx, nsply, npoints,
828 bilin, Cats, driver, mean,
829 table_name);
830 }
831 else { /* FLAG_EXT == TRUE */
832
833 /* done that earlier */
834 /*
835 int npoints_ext, *lineVect_ext = NULL;
836 double **obsVect_ext;
837 struct Point *observ_ext;
838
839 observ_ext =
840 P_Read_Vector_Region_Map(&In_ext,
841 &elaboration_reg,
842 &npoints_ext, dim_vect,
843 1);
844 */
845
846 obsVect_ext = G_alloc_matrix(npoints_ext, 3); /* Observation vector_ext */
847 lineVect_ext = G_alloc_ivector(npoints_ext);
848
849 for (i = 0; i < npoints_ext; i++) { /* Setting obsVect_ext vector & Q matrix */
850 obsVect_ext[i][0] = observ_ext[i].coordX;
851 obsVect_ext[i][1] = observ_ext[i].coordY;
852 obsVect_ext[i][2] = observ_ext[i].coordZ - mean;
853 lineVect_ext[i] = observ_ext[i].lineID;
854 }
855
856 G_free(observ_ext);
857
858 G_debug(1, "Interpolation: (%d,%d): Sparse_Points...",
859 subregion_row, subregion_col);
860 P_Sparse_Points(&Out, &elaboration_reg, general_box,
861 overlap_box, obsVect_ext, parVect,
862 lineVect_ext, stepE, stepN,
863 dims.overlap, nsplx, nsply,
864 npoints_ext, bilin, Cats, driver,
865 mean, table_name);
866
867 G_free_matrix(obsVect_ext);
868 G_free_ivector(lineVect_ext);
869 } /* END FLAG_EXT == TRUE */
870 } /* END GRID == FALSE */
871 G_free_vector(parVect);
872 G_free_matrix(obsVect);
873 G_free_ivector(lineVect);
874 }
875 else {
876 if (observ)
877 G_free(observ);
878 if (observ_ext)
879 G_free(observ_ext);
880 if (npoints == 0)
881 G_warning(_("No data within this subregion. "
882 "Consider increasing spline step values."));
883 }
884 } /*! END WHILE; last_column = TRUE */
885 } /*! END WHILE; last_row = TRUE */
886
887 G_verbose_message(_("Writing output..."));
888 /* Writing the output raster map */
889 if (grid == TRUE) {
890 int row, col;
891 DCELL *drastbuf, dval;
892
893
894 if (have_mask) {
895 Segment_close(&mask_seg); /* close segment structure */
896 }
897
898 drastbuf = Rast_allocate_buf(DCELL_TYPE);
899 for (row = 0; row < nrows; row++) {
900 G_percent(row, nrows, 2);
901 for (col = 0; col < ncols; col++) {
902 Segment_get(&out_seg, &dval, row, col);
903 drastbuf[col] = dval;
904 }
905 Rast_put_d_row(raster, drastbuf);
906 }
907
908 Rast_close(raster);
909
910 Segment_close(&out_seg); /* close segment structure */
911 /* set map title */
912 sprintf(title, "%s interpolation with Tykhonov regularization",
913 type_opt->answer);
914 Rast_put_cell_title(out_map_opt->answer, title);
915 /* write map history */
916 Rast_short_history(out_map_opt->answer, "raster", &history);
917 Rast_command_history(&history);
918 Rast_write_history(out_map_opt->answer, &history);
919 }
920 /* Writing to the output vector map the points from the overlapping zones */
921 else if (flag_auxiliar == TRUE) {
922 if (ext == FALSE)
923 P_Aux_to_Vector(&In, &Out, driver, table_name);
924 else
925 P_Aux_to_Vector(&In_ext, &Out, driver, table_name);
926
927 /* Drop auxiliary table */
928 G_debug(1, "%s: Dropping <%s>", argv[0], table_name);
929 if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
930 G_fatal_error(_("Auxiliary table could not be dropped"));
931 }
932
933 db_close_database_shutdown_driver(driver);
934
935 Vect_close(&In);
936 if (ext != FALSE)
937 Vect_close(&In_ext);
938 if (vector)
939 Vect_close(&Out);
940
941 G_done_msg(" ");
942
943 exit(EXIT_SUCCESS);
944} /*END MAIN */
Note: See TracBrowser for help on using the repository browser.