← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to vector/v.outlier/main.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
vector/v.outlier/main.c
 
1
 
 
2
/**********************************************************************
 
3
 *
 
4
 * MODULE:       v.outlier
 
5
 *
 
6
 * AUTHOR(S):    Roberto Antolin
 
7
 *               
 
8
 * PURPOSE:      Removal of data outliers
 
9
 *
 
10
 * COPYRIGHT:    (C) 2006 by Politecnico di Milano -
 
11
 *                           Polo Regionale di Como
 
12
 *
 
13
 *               This program is free software under the
 
14
 *               GNU General Public License (>=v2).
 
15
 *               Read the file COPYING that comes with GRASS
 
16
 *               for details.
 
17
 *
 
18
 **********************************************************************/
 
19
 
 
20
/* INCLUDES */
 
21
#include <stdlib.h>
 
22
#include <string.h>
 
23
#include <math.h>
 
24
#include "outlier.h"
 
25
 
 
26
/* GLOBAL VARIABLES */
 
27
int nsply, nsplx;
 
28
double stepN, stepE, Thres_Outlier;
 
29
 
 
30
/*--------------------------------------------------------------------*/
 
31
int main(int argc, char *argv[])
 
32
{
 
33
    /* Variables declarations */
 
34
    int nsplx_adj, nsply_adj;
 
35
    int nsubregion_col, nsubregion_row;
 
36
    int subregion = 0, nsubregions = 0;
 
37
    double N_extension, E_extension, edgeE, edgeN;
 
38
    int dim_vect, nparameters, BW, npoints;
 
39
    double mean, lambda;
 
40
    const char *dvr, *db, *mapset;
 
41
    char table_name[GNAME_MAX];
 
42
    char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
 
43
 
 
44
    int last_row, last_column, flag_auxiliar = FALSE;
 
45
    int filter_mode;
 
46
 
 
47
    int *lineVect;
 
48
    double *TN, *Q, *parVect;   /* Interpolating and least-square vectors */
 
49
    double **N, **obsVect;      /* Interpolation and least-square matrix */
 
50
 
 
51
    /* Structs declarations */
 
52
    struct Map_info In, Out, Outlier, Qgis;
 
53
    struct Option *in_opt, *out_opt, *outlier_opt, *qgis_opt, *stepE_opt,
 
54
        *stepN_opt, *lambda_f_opt, *Thres_O_opt, *filter_opt;
 
55
    struct Flag *spline_step_flag;
 
56
    struct GModule *module;
 
57
 
 
58
    struct Reg_dimens dims;
 
59
    struct Cell_head elaboration_reg, original_reg;
 
60
    struct bound_box general_box, overlap_box;
 
61
 
 
62
    struct Point *observ;
 
63
 
 
64
    dbDriver *driver;
 
65
 
 
66
    /*----------------------------------------------------------------*/
 
67
    /* Options declaration */
 
68
    module = G_define_module();
 
69
    G_add_keyword(_("vector"));
 
70
    G_add_keyword(_("statistics"));
 
71
    module->description = _("Removes outliers from vector point data.");
 
72
 
 
73
    spline_step_flag = G_define_flag();
 
74
    spline_step_flag->key = 'e';
 
75
    spline_step_flag->label = _("Estimate point density and distance");
 
76
    spline_step_flag->description =
 
77
        _("Estimate point density and distance for the input vector points within the current region extends and quit");
 
78
 
 
79
    in_opt = G_define_standard_option(G_OPT_V_INPUT);
 
80
 
 
81
    out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
 
82
 
 
83
    outlier_opt = G_define_option();
 
84
    outlier_opt->key = "outlier";
 
85
    outlier_opt->type = TYPE_STRING;
 
86
    outlier_opt->key_desc = "name";
 
87
    outlier_opt->required = YES;
 
88
    outlier_opt->gisprompt = "new,vector,vector";
 
89
    outlier_opt->description = _("Name of output outlier vector map");
 
90
 
 
91
    qgis_opt = G_define_option();
 
92
    qgis_opt->key = "qgis";
 
93
    qgis_opt->type = TYPE_STRING;
 
94
    qgis_opt->key_desc = "name";
 
95
    qgis_opt->required = NO;
 
96
    qgis_opt->gisprompt = "new,vector,vector";
 
97
    qgis_opt->description = _("Name of vector map for visualization in QGIS");
 
98
 
 
99
    stepE_opt = G_define_option();
 
100
    stepE_opt->key = "ew_step";
 
101
    stepE_opt->type = TYPE_DOUBLE;
 
102
    stepE_opt->required = NO;
 
103
    stepE_opt->answer = "10";
 
104
    stepE_opt->description =
 
105
        _("Length of each spline step in the east-west direction");
 
106
    stepE_opt->guisection = _("Settings");
 
107
 
 
108
    stepN_opt = G_define_option();
 
109
    stepN_opt->key = "ns_step";
 
110
    stepN_opt->type = TYPE_DOUBLE;
 
111
    stepN_opt->required = NO;
 
112
    stepN_opt->answer = "10";
 
113
    stepN_opt->description =
 
114
        _("Length of each spline step in the north-south direction");
 
115
    stepN_opt->guisection = _("Settings");
 
116
 
 
117
    lambda_f_opt = G_define_option();
 
118
    lambda_f_opt->key = "lambda";
 
119
    lambda_f_opt->type = TYPE_DOUBLE;
 
120
    lambda_f_opt->required = NO;
 
121
    lambda_f_opt->description = _("Tykhonov regularization weight");
 
122
    lambda_f_opt->answer = "0.1";
 
123
    lambda_f_opt->guisection = _("Settings");
 
124
 
 
125
    Thres_O_opt = G_define_option();
 
126
    Thres_O_opt->key = "threshold";
 
127
    Thres_O_opt->type = TYPE_DOUBLE;
 
128
    Thres_O_opt->required = NO;
 
129
    Thres_O_opt->description = _("Threshold for the outliers");
 
130
    Thres_O_opt->answer = "50";
 
131
 
 
132
    filter_opt = G_define_option();
 
133
    filter_opt->key = "filter";
 
134
    filter_opt->type = TYPE_STRING;
 
135
    filter_opt->required = NO;
 
136
    filter_opt->description = _("Filtering option");
 
137
    filter_opt->options = "both,positive,negative";
 
138
    filter_opt->answer = "both";
 
139
 
 
140
    /* Parsing */
 
141
    G_gisinit(argv[0]);
 
142
    if (G_parser(argc, argv))
 
143
        exit(EXIT_FAILURE);
 
144
 
 
145
    if (!(db = G_getenv_nofatal2("DB_DATABASE", G_VAR_MAPSET)))
 
146
        G_fatal_error(_("Unable to read name of database"));
 
147
 
 
148
    if (!(dvr = G_getenv_nofatal2("DB_DRIVER", G_VAR_MAPSET)))
 
149
        G_fatal_error(_("Unable to read name of driver"));
 
150
 
 
151
    stepN = atof(stepN_opt->answer);
 
152
    stepE = atof(stepE_opt->answer);
 
153
    lambda = atof(lambda_f_opt->answer);
 
154
    Thres_Outlier = atof(Thres_O_opt->answer);
 
155
 
 
156
    filter_mode = 0;
 
157
    if (strcmp(filter_opt->answer, "positive") == 0)
 
158
        filter_mode = 1;
 
159
    else if (strcmp(filter_opt->answer, "negative") == 0)
 
160
        filter_mode = -1;
 
161
    P_set_outlier_fn(filter_mode);
 
162
 
 
163
    flag_auxiliar = FALSE;
 
164
 
 
165
    /* Checking vector names */
 
166
    Vect_check_input_output_name(in_opt->answer, out_opt->answer,
 
167
                                 G_FATAL_EXIT);
 
168
 
 
169
    if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL) {
 
170
        G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
 
171
    }
 
172
 
 
173
    /* Setting auxiliar table's name */
 
174
    if (G_name_is_fully_qualified(out_opt->answer, xname, xmapset)) {
 
175
        sprintf(table_name, "%s_aux", xname);
 
176
    }
 
177
    else
 
178
        sprintf(table_name, "%s_aux", out_opt->answer);
 
179
 
 
180
    /* Something went wrong in a previous v.outlier execution */
 
181
    if (db_table_exists(dvr, db, table_name)) {
 
182
        /* Start driver and open db */
 
183
        driver = db_start_driver_open_database(dvr, db);
 
184
        if (driver == NULL)
 
185
            G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
 
186
                          dvr);
 
187
        db_set_error_handler_driver(driver);
 
188
 
 
189
        if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
 
190
            G_fatal_error(_("Old auxiliar table could not be dropped"));
 
191
        db_close_database_shutdown_driver(driver);
 
192
    }
 
193
 
 
194
    /* Open input vector */
 
195
    Vect_set_open_level(1);     /* WITHOUT TOPOLOGY */
 
196
    if (1 > Vect_open_old(&In, in_opt->answer, mapset))
 
197
        G_fatal_error(_("Unable to open vector map <%s> at the topological level"),
 
198
                      in_opt->answer);
 
199
 
 
200
    /* Input vector must be 3D */
 
201
    if (!Vect_is_3d(&In))
 
202
        G_fatal_error(_("Input vector map <%s> is not 3D!"), in_opt->answer);
 
203
 
 
204
    /* Estimate point density and mean distance for current region */
 
205
    if (spline_step_flag->answer) {
 
206
        double dens, dist;
 
207
        if (P_estimate_splinestep(&In, &dens, &dist) == 0) {
 
208
            G_message("Estimated point density: %.4g", dens);
 
209
            G_message("Estimated mean distance between points: %.4g", dist);
 
210
        }
 
211
        else
 
212
            G_warning(_("No points in current region!"));
 
213
        
 
214
        Vect_close(&In);
 
215
        exit(EXIT_SUCCESS);
 
216
    }
 
217
 
 
218
    /* Open output vector */
 
219
    if (qgis_opt->answer)
 
220
        if (0 > Vect_open_new(&Qgis, qgis_opt->answer, WITHOUT_Z))
 
221
            G_fatal_error(_("Unable to create vector map <%s>"),
 
222
                          qgis_opt->answer);
 
223
 
 
224
    if (0 > Vect_open_new(&Out, out_opt->answer, WITH_Z)) {
 
225
        Vect_close(&Qgis);
 
226
        G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
 
227
    }
 
228
 
 
229
    if (0 > Vect_open_new(&Outlier, outlier_opt->answer, WITH_Z)) {
 
230
        Vect_close(&Out);
 
231
        Vect_close(&Qgis);
 
232
        G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
 
233
    }
 
234
 
 
235
    /* Copy vector Head File */
 
236
    Vect_copy_head_data(&In, &Out);
 
237
    Vect_hist_copy(&In, &Out);
 
238
    Vect_hist_command(&Out);
 
239
 
 
240
    Vect_copy_head_data(&In, &Outlier);
 
241
    Vect_hist_copy(&In, &Outlier);
 
242
    Vect_hist_command(&Outlier);
 
243
 
 
244
    if (qgis_opt->answer) {
 
245
        Vect_copy_head_data(&In, &Qgis);
 
246
        Vect_hist_copy(&In, &Qgis);
 
247
        Vect_hist_command(&Qgis);
 
248
    }
 
249
 
 
250
    /* Open driver and database */
 
251
    driver = db_start_driver_open_database(dvr, db);
 
252
    if (driver == NULL)
 
253
        G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
 
254
                      dvr);
 
255
    db_set_error_handler_driver(driver);
 
256
 
 
257
    /* Create auxiliar table */
 
258
    if ((flag_auxiliar =
 
259
         P_Create_Aux2_Table(driver, table_name)) == FALSE)
 
260
        G_fatal_error(_("It was impossible to create <%s> table."), table_name);
 
261
 
 
262
    db_create_index2(driver, table_name, "ID");
 
263
    /* sqlite likes that ??? */
 
264
    db_close_database_shutdown_driver(driver);
 
265
    driver = db_start_driver_open_database(dvr, db);
 
266
 
 
267
    /* Setting regions and boxes */
 
268
    G_get_set_window(&original_reg);
 
269
    G_get_set_window(&elaboration_reg);
 
270
    Vect_region_box(&elaboration_reg, &overlap_box);
 
271
    Vect_region_box(&elaboration_reg, &general_box);
 
272
 
 
273
    /*------------------------------------------------------------------
 
274
      | Subdividing and working with tiles:                                                                     
 
275
      | Each original region will be divided into several subregions. 
 
276
      | Each one will be overlaped by its neighbouring subregions. 
 
277
      | The overlapping is calculated as a fixed OVERLAP_SIZE times
 
278
      | the largest spline step plus 2 * edge
 
279
      ----------------------------------------------------------------*/
 
280
 
 
281
    /* Fixing parameters of the elaboration region */
 
282
    P_zero_dim(&dims);          /* Set dim struct to zero */
 
283
 
 
284
    nsplx_adj = NSPLX_MAX;
 
285
    nsply_adj = NSPLY_MAX;
 
286
    if (stepN > stepE)
 
287
        dims.overlap = OVERLAP_SIZE * stepN;
 
288
    else
 
289
        dims.overlap = OVERLAP_SIZE * stepE;
 
290
    P_get_edge(P_BILINEAR, &dims, stepE, stepN);
 
291
    P_set_dim(&dims, stepE, stepN, &nsplx_adj, &nsply_adj);
 
292
 
 
293
    G_verbose_message(_("Adjusted EW splines %d"), nsplx_adj);
 
294
    G_verbose_message(_("Adjusted NS splines %d"), nsply_adj);
 
295
 
 
296
    /* calculate number of subregions */
 
297
    edgeE = dims.ew_size - dims.overlap - 2 * dims.edge_v;
 
298
    edgeN = dims.sn_size - dims.overlap - 2 * dims.edge_h;
 
299
 
 
300
    N_extension = original_reg.north - original_reg.south;
 
301
    E_extension = original_reg.east - original_reg.west;
 
302
 
 
303
    nsubregion_col = ceil(E_extension / edgeE) + 0.5;
 
304
    nsubregion_row = ceil(N_extension / edgeN) + 0.5;
 
305
 
 
306
    if (nsubregion_col < 0)
 
307
        nsubregion_col = 0;
 
308
    if (nsubregion_row < 0)
 
309
        nsubregion_row = 0;
 
310
 
 
311
    nsubregions = nsubregion_row * nsubregion_col;
 
312
 
 
313
    elaboration_reg.south = original_reg.north;
 
314
    last_row = FALSE;
 
315
 
 
316
    while (last_row == FALSE) { /* For each row */
 
317
 
 
318
        P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 
319
                      GENERAL_ROW);
 
320
 
 
321
        if (elaboration_reg.north > original_reg.north) {       /* First row */
 
322
 
 
323
            P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 
324
                          FIRST_ROW);
 
325
        }
 
326
 
 
327
        if (elaboration_reg.south <= original_reg.south) {      /* Last row */
 
328
 
 
329
            P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 
330
                          LAST_ROW);
 
331
            last_row = TRUE;
 
332
        }
 
333
 
 
334
        nsply =
 
335
            ceil((elaboration_reg.north -
 
336
                  elaboration_reg.south) / stepN) + 0.5;
 
337
        /*
 
338
        if (nsply > NSPLY_MAX)
 
339
            nsply = NSPLY_MAX;
 
340
        */
 
341
        G_debug(1, "nsply = %d", nsply);
 
342
 
 
343
        elaboration_reg.east = original_reg.west;
 
344
        last_column = FALSE;
 
345
 
 
346
        while (last_column == FALSE) {  /* For each column */
 
347
 
 
348
            subregion++;
 
349
            if (nsubregions > 1)
 
350
                G_message(_("Processing subregion %d of %d..."), subregion, nsubregions);
 
351
            else /* v.outlier -e will report mean point distance: */
 
352
                G_warning(_("No subregions found! Check values for 'ew_step' and 'ns_step' parameters"));
 
353
 
 
354
            P_set_regions(&elaboration_reg, &general_box, &overlap_box, dims,
 
355
                          GENERAL_COLUMN);
 
356
 
 
357
            if (elaboration_reg.west < original_reg.west) {     /* First column */
 
358
 
 
359
                P_set_regions(&elaboration_reg, &general_box, &overlap_box,
 
360
                              dims, FIRST_COLUMN);
 
361
            }
 
362
 
 
363
            if (elaboration_reg.east >= original_reg.east) {    /* Last column */
 
364
 
 
365
                P_set_regions(&elaboration_reg, &general_box, &overlap_box,
 
366
                              dims, LAST_COLUMN);
 
367
                last_column = TRUE;
 
368
            }
 
369
            nsplx =
 
370
                ceil((elaboration_reg.east -
 
371
                      elaboration_reg.west) / stepE) + 0.5;
 
372
            /*
 
373
            if (nsplx > NSPLX_MAX)
 
374
                nsplx = NSPLX_MAX;
 
375
            */
 
376
            G_debug(1, "nsplx = %d", nsplx);
 
377
 
 
378
            /*Setting the active region */
 
379
            dim_vect = nsplx * nsply;
 
380
            observ =
 
381
                P_Read_Vector_Region_Map(&In, &elaboration_reg, &npoints,
 
382
                                         dim_vect, 1);
 
383
 
 
384
            if (npoints > 0) {  /* If there is any point falling into elaboration_reg... */
 
385
                int i;
 
386
 
 
387
                nparameters = nsplx * nsply;
 
388
 
 
389
                /* Mean calculation */
 
390
                mean = P_Mean_Calc(&elaboration_reg, observ, npoints);
 
391
 
 
392
                /* Least Squares system */
 
393
                G_debug(1, "Allocation memory for bilinear interpolation");
 
394
                BW = P_get_BandWidth(P_BILINEAR, nsply);        /* Bilinear interpolation */
 
395
                N = G_alloc_matrix(nparameters, BW);    /* Normal matrix */
 
396
                TN = G_alloc_vector(nparameters);       /* vector */
 
397
                parVect = G_alloc_vector(nparameters);  /* Bicubic parameters vector */
 
398
                obsVect = G_alloc_matrix(npoints, 3);   /* Observation vector */
 
399
                Q = G_alloc_vector(npoints);    /* "a priori" var-cov matrix */
 
400
                lineVect = G_alloc_ivector(npoints);
 
401
 
 
402
                /* Setting obsVect vector & Q matrix */
 
403
                for (i = 0; i < npoints; i++) {
 
404
                    obsVect[i][0] = observ[i].coordX;
 
405
                    obsVect[i][1] = observ[i].coordY;
 
406
                    obsVect[i][2] = observ[i].coordZ - mean;
 
407
                    lineVect[i] = observ[i].lineID;
 
408
                    Q[i] = 1;   /* Q=I */
 
409
                }
 
410
 
 
411
                G_free(observ);
 
412
 
 
413
                G_verbose_message(_("Bilinear interpolation"));
 
414
                normalDefBilin(N, TN, Q, obsVect, stepE, stepN, nsplx,
 
415
                               nsply, elaboration_reg.west,
 
416
                               elaboration_reg.south, npoints, nparameters,
 
417
                               BW);
 
418
                nCorrectGrad(N, lambda, nsplx, nsply, stepE, stepN);
 
419
                G_math_solver_cholesky_sband(N, parVect, TN, nparameters, BW);
 
420
 
 
421
                G_free_matrix(N);
 
422
                G_free_vector(TN);
 
423
                G_free_vector(Q);
 
424
 
 
425
                G_verbose_message(_("Outlier detection"));
 
426
                if (qgis_opt->answer)
 
427
                    P_Outlier(&Out, &Outlier, &Qgis, elaboration_reg,
 
428
                              general_box, overlap_box, obsVect, parVect,
 
429
                              mean, dims.overlap, lineVect, npoints,
 
430
                              driver, table_name);
 
431
                else
 
432
                    P_Outlier(&Out, &Outlier, NULL, elaboration_reg,
 
433
                              general_box, overlap_box, obsVect, parVect,
 
434
                              mean, dims.overlap, lineVect, npoints,
 
435
                              driver, table_name);
 
436
 
 
437
 
 
438
                G_free_vector(parVect);
 
439
                G_free_matrix(obsVect);
 
440
                G_free_ivector(lineVect);
 
441
 
 
442
            }                   /*! END IF; npoints > 0 */
 
443
            else {
 
444
                G_free(observ);
 
445
                G_warning(_("No data within this subregion. "
 
446
                            "Consider increasing spline step values."));
 
447
            }
 
448
        }                       /*! END WHILE; last_column = TRUE */
 
449
    }                           /*! END WHILE; last_row = TRUE */
 
450
 
 
451
    /* Drop auxiliar table */
 
452
    if (npoints > 0) {
 
453
        G_debug(1, "%s: Dropping <%s>", argv[0], table_name);
 
454
        if (P_Drop_Aux_Table(driver, table_name) != DB_OK)
 
455
            G_fatal_error(_("Auxiliary table could not be dropped"));
 
456
    }
 
457
 
 
458
    db_close_database_shutdown_driver(driver);
 
459
 
 
460
    Vect_close(&In);
 
461
    Vect_close(&Out);
 
462
    Vect_close(&Outlier);
 
463
    if (qgis_opt->answer) {
 
464
        Vect_build(&Qgis);
 
465
        Vect_close(&Qgis);
 
466
    }
 
467
 
 
468
    G_done_msg(" ");
 
469
 
 
470
    exit(EXIT_SUCCESS);
 
471
}                               /*END MAIN */