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

« back to all changes in this revision

Viewing changes to vector/v.lidar.growing/main.c

Tags: 7.0.0~rc1+ds1-1~exp1
* 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:
 
1
 
 
2
/****************************************************************
 
3
 *                                                              *
 
4
 * MODULE:      v.lidar.growing                         *
 
5
 *                                                              *
 
6
 * AUTHOR(S):   Roberto Antolin & Gonzalo Moreno                *
 
7
 *                                                              *
 
8
 * PURPOSE:     Building contour determination and Region       * 
 
9
 *              Growing algorithm for determining the building  *
 
10
 *              inside                                          *
 
11
 *                                                              *
 
12
 * COPYRIGHT:   (C) 2006 by Politecnico di Milano -             *
 
13
 *                           Polo Regionale di Como             *
 
14
 *                                                              *
 
15
 *               This program is free software under the        *
 
16
 *               GNU General Public License (>=v2).             *
 
17
 *               Read the file COPYING that comes with GRASS    *
 
18
 *               for details.                                   *
 
19
 *                                                              *
 
20
 ****************************************************************/
 
21
 
 
22
 /*INCLUDES*/
 
23
#include <stdlib.h>
 
24
#include <string.h>
 
25
#include <math.h>
 
26
#include "growing.h"
 
27
 
 
28
    /* GLOBAL DEFINITIONS */
 
29
int nsply, nsplx, count_obj;
 
30
double stepN, stepE;
 
31
 
 
32
/*--------------------------------------------------------------------------------*/
 
33
int main(int argc, char *argv[])
 
34
{
 
35
 
 
36
    /* Variables' declarations */
 
37
    int row, nrows, col, ncols, MaxPoints;
 
38
    int nsubregion_col, nsubregion_row;
 
39
    int subregion = 0, nsubregions = 0;
 
40
    int last_row, last_column;
 
41
    int nlines, nlines_first, line_num;
 
42
    int more;
 
43
    int clas, region = TRUE;
 
44
    double Z_interp;
 
45
    double Thres_j, Thres_d, ew_resol, ns_resol;
 
46
    double minNS, minEW, maxNS, maxEW;
 
47
    const char *mapset;
 
48
    char buf[1024], table_name[GNAME_MAX];
 
49
    char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
 
50
 
 
51
    int colorBordo, ripieno, conta, lungPunti, lungHull, xi, c1, c2;
 
52
    double altPiano;
 
53
    extern double **P, **cvxHull, **punti_bordo;
 
54
 
 
55
    /* Struct declarations */
 
56
    struct Cell_head elaboration_reg, original_reg;
 
57
    struct element_grow **raster_matrix;
 
58
 
 
59
    struct Map_info In, Out, First;
 
60
    struct Option *in_opt, *out_opt, *first_opt, *Thres_j_opt, *Thres_d_opt;
 
61
    struct GModule *module;
 
62
 
 
63
    struct line_pnts *points, *points_first;
 
64
    struct line_cats *Cats, *Cats_first;
 
65
 
 
66
    struct field_info *field;
 
67
    dbDriver *driver;
 
68
    dbString sql;
 
69
    dbTable *table;
 
70
    dbCursor cursor;
 
71
 
 
72
/*------------------------------------------------------------------------------------------*/
 
73
    /* Options' declaration */ ;
 
74
    module = G_define_module();
 
75
    G_add_keyword(_("vector"));
 
76
    G_add_keyword(_("LIDAR"));
 
77
    module->description =
 
78
        _("Building contour determination and Region Growing "
 
79
          "algorithm for determining the building inside");
 
80
 
 
81
    in_opt = G_define_standard_option(G_OPT_V_INPUT);
 
82
    in_opt->description =
 
83
        _("Input vector (v.lidar.edgedetection output)");
 
84
 
 
85
    out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
 
86
 
 
87
    first_opt = G_define_option();
 
88
    first_opt->key = "first";
 
89
    first_opt->type = TYPE_STRING;
 
90
    first_opt->key_desc = "name";
 
91
    first_opt->required = YES;
 
92
    first_opt->gisprompt = "old,vector,vector";
 
93
    first_opt->description = _("Name of the first pulse vector map");
 
94
 
 
95
    Thres_j_opt = G_define_option();
 
96
    Thres_j_opt->key = "tj";
 
97
    Thres_j_opt->type = TYPE_DOUBLE;
 
98
    Thres_j_opt->required = NO;
 
99
    Thres_j_opt->description =
 
100
        _("Threshold for cell object frequency in region growing");
 
101
    Thres_j_opt->answer = "0.2";
 
102
 
 
103
    Thres_d_opt = G_define_option();
 
104
    Thres_d_opt->key = "td";
 
105
    Thres_d_opt->type = TYPE_DOUBLE;
 
106
    Thres_d_opt->required = NO;
 
107
    Thres_d_opt->description =
 
108
        _("Threshold for double pulse in region growing");
 
109
    Thres_d_opt->answer = "0.6";
 
110
 
 
111
    /* Parsing */
 
112
    G_gisinit(argv[0]);
 
113
    if (G_parser(argc, argv))
 
114
        exit(EXIT_FAILURE);
 
115
 
 
116
    Thres_j = atof(Thres_j_opt->answer);
 
117
    Thres_d = atof(Thres_d_opt->answer);
 
118
 
 
119
    Thres_j += 1;
 
120
 
 
121
    /* Open input vector */
 
122
    Vect_check_input_output_name(in_opt->answer, out_opt->answer,
 
123
                                 G_FATAL_EXIT);
 
124
 
 
125
    if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL) {
 
126
        G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
 
127
    }
 
128
 
 
129
    /* Setting auxiliar table's name */
 
130
    if (G_name_is_fully_qualified(in_opt->answer, xname, xmapset)) {
 
131
        sprintf(table_name, "%s_edge_Interpolation", xname);
 
132
    }
 
133
    else
 
134
        sprintf(table_name, "%s_edge_Interpolation", in_opt->answer);
 
135
 
 
136
    Vect_set_open_level(1);     /* WITHOUT TOPOLOGY */
 
137
    if (Vect_open_old(&In, in_opt->answer, mapset) < 1)
 
138
        G_fatal_error(_("Unable to open vector map <%s>"), in_opt->answer);
 
139
 
 
140
    Vect_set_open_level(1);     /* WITHOUT TOPOLOGY */
 
141
    if (Vect_open_old(&First, first_opt->answer, mapset) < 1)
 
142
        G_fatal_error(_("Unable to open vector map <%s>"), first_opt->answer);
 
143
 
 
144
    /* Open output vector */
 
145
    if (0 > Vect_open_new(&Out, out_opt->answer, WITH_Z)) {
 
146
        Vect_close(&In);
 
147
        Vect_close(&First);
 
148
        exit(EXIT_FAILURE);
 
149
    }
 
150
 
 
151
    /* Copy vector Head File */
 
152
    Vect_copy_head_data(&In, &Out);
 
153
    Vect_hist_copy(&In, &Out);
 
154
    Vect_hist_command(&Out);
 
155
 
 
156
    /* Starting driver and open db for edgedetection interpolation table */
 
157
    field = Vect_get_field(&In, F_INTERPOLATION);
 
158
    /*if (field == NULL)
 
159
       G_fatal_error (_("Cannot read field info")); */
 
160
 
 
161
    driver = db_start_driver_open_database(field->driver, field->database);
 
162
    if (driver == NULL)
 
163
        G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
 
164
                      field->driver);
 
165
 
 
166
    /* is this the right place to open the cursor ??? */
 
167
    
 
168
    db_init_string(&sql);
 
169
    db_zero_string(&sql);
 
170
 
 
171
    sprintf(buf, "SELECT Interp,ID FROM %s", table_name);
 
172
    G_debug(1, "buf: %s", buf);
 
173
    db_append_string(&sql, buf);
 
174
 
 
175
    if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
 
176
        G_fatal_error(_("Unable to open table <%s>"), table_name);
 
177
 
 
178
    count_obj = 1;
 
179
 
 
180
    /* no topology, get number of lines in input vector */
 
181
    nlines = 0;
 
182
    points = Vect_new_line_struct();
 
183
    Cats = Vect_new_cats_struct();
 
184
    Vect_rewind(&In);
 
185
    while (Vect_read_next_line(&In, points, Cats) > 0) {
 
186
        nlines++;
 
187
    }
 
188
    Vect_rewind(&In);
 
189
 
 
190
    /* no topology, get number of lines in first pulse input vector */
 
191
    nlines_first = 0;
 
192
    points_first = Vect_new_line_struct();
 
193
    Cats_first = Vect_new_cats_struct();
 
194
    Vect_rewind(&First);
 
195
    while (Vect_read_next_line(&First, points_first, Cats_first) > 0) {
 
196
        nlines_first++;
 
197
    }
 
198
    Vect_rewind(&First);
 
199
 
 
200
    /* Setting regions and boxes */
 
201
    G_debug(1, _("Setting regions and boxes"));
 
202
    G_get_set_window(&original_reg);
 
203
    G_get_set_window(&elaboration_reg);
 
204
 
 
205
    /*  Fixing parameters of the elaboration region */
 
206
    /*! The original_region will be divided into subregions */
 
207
    ew_resol = original_reg.ew_res;
 
208
    ns_resol = original_reg.ns_res;
 
209
 
 
210
    /* calculate number of subregions */
 
211
    nsubregion_col = ceil((original_reg.east - original_reg.west) / (LATO * ew_resol)) + 0.5;
 
212
    nsubregion_row = ceil((original_reg.north - original_reg.south) / (LATO * ns_resol)) + 0.5;
 
213
 
 
214
    if (nsubregion_col < 0)
 
215
        nsubregion_col = 0;
 
216
    if (nsubregion_row < 0)
 
217
        nsubregion_row = 0;
 
218
 
 
219
    nsubregions = nsubregion_row * nsubregion_col;
 
220
 
 
221
    /* Subdividing and working with tiles */
 
222
    elaboration_reg.south = original_reg.north;
 
223
    last_row = FALSE;
 
224
 
 
225
    while (last_row == FALSE) { /* For each strip of LATO rows */
 
226
 
 
227
        elaboration_reg.north = elaboration_reg.south;
 
228
 
 
229
        if (elaboration_reg.north > original_reg.north) /* First row */
 
230
            elaboration_reg.north = original_reg.north;
 
231
 
 
232
        elaboration_reg.south = elaboration_reg.north - LATO * ns_resol;
 
233
        if (elaboration_reg.south <= original_reg.south) {      /* Last row */
 
234
            elaboration_reg.south = original_reg.south;
 
235
            last_row = TRUE;
 
236
        }
 
237
 
 
238
        elaboration_reg.east = original_reg.west;
 
239
        last_column = FALSE;
 
240
 
 
241
        while (last_column == FALSE) {  /* For each strip of LATO columns */
 
242
            struct bound_box elaboration_box;
 
243
 
 
244
            subregion++;
 
245
            if (nsubregions > 1)
 
246
                G_message(_("subregion %d of %d"), subregion, nsubregions);
 
247
 
 
248
            elaboration_reg.west = elaboration_reg.east;
 
249
            if (elaboration_reg.west < original_reg.west)       /* First column */
 
250
                elaboration_reg.west = original_reg.west;
 
251
 
 
252
            elaboration_reg.east = elaboration_reg.west + LATO * ew_resol;
 
253
 
 
254
            if (elaboration_reg.east >= original_reg.east) {    /* Last column */
 
255
                elaboration_reg.east = original_reg.east;
 
256
                last_column = TRUE;
 
257
            }
 
258
 
 
259
            /* Setting the active region */
 
260
            elaboration_reg.ns_res = ns_resol;
 
261
            elaboration_reg.ew_res = ew_resol;
 
262
            nrows = (elaboration_reg.north - elaboration_reg.south) / ns_resol + 0.1;
 
263
            ncols = (elaboration_reg.east - elaboration_reg.west) / ew_resol + 0.1;
 
264
            elaboration_reg.rows = nrows;
 
265
            elaboration_reg.cols = ncols;
 
266
 
 
267
            G_debug(1, _("Rows = %d"), nrows);
 
268
            G_debug(1, _("Columns = %d"), ncols);
 
269
 
 
270
            raster_matrix = structMatrix(0, nrows, 0, ncols);
 
271
            MaxPoints = nrows * ncols;
 
272
 
 
273
            /* Initializing matrix */
 
274
            for (row = 0; row <= nrows; row++) {
 
275
                for (col = 0; col <= ncols; col++) {
 
276
                    raster_matrix[row][col].interp = 0;
 
277
                    raster_matrix[row][col].fi = 0;
 
278
                    raster_matrix[row][col].bordo = 0;
 
279
                    raster_matrix[row][col].dueImp = SINGLE_PULSE;
 
280
                    raster_matrix[row][col].orig = 0;
 
281
                    raster_matrix[row][col].fo = 0;
 
282
                    raster_matrix[row][col].clas = PRE_TERRAIN;
 
283
                    raster_matrix[row][col].fc = 0;
 
284
                    raster_matrix[row][col].obj = 0;
 
285
                }
 
286
            }
 
287
 
 
288
            G_verbose_message(_("read points in input vector"));
 
289
            Vect_region_box(&elaboration_reg, &elaboration_box);
 
290
            line_num = 0;
 
291
            Vect_rewind(&In);
 
292
            while (Vect_read_next_line(&In, points, Cats) > 0) {
 
293
                line_num++;
 
294
 
 
295
                if ((Vect_point_in_box
 
296
                     (points->x[0], points->y[0], points->z[0],
 
297
                      &elaboration_box)) &&
 
298
                    ((points->x[0] != elaboration_reg.west) ||
 
299
                     (points->x[0] == original_reg.west)) &&
 
300
                    ((points->y[0] != elaboration_reg.north) ||
 
301
                     (points->y[0] == original_reg.north))) {
 
302
 
 
303
                    row =
 
304
                        (int)(Rast_northing_to_row
 
305
                              (points->y[0], &elaboration_reg));
 
306
                    col =
 
307
                        (int)(Rast_easting_to_col
 
308
                              (points->x[0], &elaboration_reg));
 
309
 
 
310
                    Z_interp = 0;
 
311
                    /* TODO: make sure the current db_fetch() usage works */
 
312
                    /* why not: */
 
313
                    /*
 
314
                    db_init_string(&sql);
 
315
                    sprintf(buf, "SELECT Interp,ID FROM %s WHERE ID=%d", table_name, line_num);
 
316
                    db_append_string(&sql, buf);
 
317
 
 
318
                    if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
 
319
                        G_fatal_error(_("Unable to open table <%s>"), table_name);
 
320
 
 
321
                    while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
 
322
                        dbColumn *Z_Interp_col;
 
323
                        dbValue *Z_Interp_value;
 
324
                        table = db_get_cursor_table(&cursor);
 
325
 
 
326
                        Z_Interp_col = db_get_table_column(table, 1);
 
327
 
 
328
                        if (db_sqltype_to_Ctype(db_get_column_sqltype(Z_Interp_col)) ==
 
329
                            DB_C_TYPE_DOUBLE)
 
330
                            Z_Interp_value = db_get_column_value(Z_Interp_col);
 
331
                        else
 
332
                            continue;
 
333
 
 
334
                        Z_interp = db_get_value_double(Z_Interp_value);
 
335
                        break;
 
336
                    }
 
337
                    db_close_cursor(&cursor);
 
338
                    db_free_string(&sql);
 
339
                    */
 
340
                    /* instead of */
 
341
                    while (1) {
 
342
                        if (db_fetch(&cursor, DB_NEXT, &more) != DB_OK ||
 
343
                            !more)
 
344
                            break;
 
345
                        dbColumn *Z_Interp_col, *ID_col;
 
346
                        dbValue *Z_Interp_value, *ID_value;
 
347
 
 
348
                        table = db_get_cursor_table(&cursor);
 
349
 
 
350
                        ID_col = db_get_table_column(table, 1);
 
351
                        if (db_sqltype_to_Ctype(db_get_column_sqltype(ID_col))
 
352
                            == DB_C_TYPE_INT)
 
353
                            ID_value = db_get_column_value(ID_col);
 
354
                        else
 
355
                            continue;
 
356
 
 
357
                        if (db_get_value_int(ID_value) == line_num) {
 
358
                            Z_Interp_col = db_get_table_column(table, 0);
 
359
                            if (db_sqltype_to_Ctype
 
360
                                (db_get_column_sqltype(Z_Interp_col)) ==
 
361
                                DB_C_TYPE_DOUBLE)
 
362
                                Z_Interp_value =
 
363
                                    db_get_column_value(Z_Interp_col);
 
364
                            else
 
365
                                continue;
 
366
                            Z_interp = db_get_value_double(Z_Interp_value);
 
367
                            break;
 
368
                        }
 
369
                    }
 
370
 
 
371
                    raster_matrix[row][col].interp += Z_interp;
 
372
                    raster_matrix[row][col].fi++;
 
373
 
 
374
                    /*if (( clas = Vect_get_line_cat (&In, line_num, F_EDGE_DETECTION_CLASS) ) != UNKNOWN_EDGE) { */
 
375
                    if (Vect_cat_get(Cats, F_EDGE_DETECTION_CLASS, &clas)) {
 
376
                        raster_matrix[row][col].clas += clas;
 
377
                        raster_matrix[row][col].fc++;
 
378
                    }
 
379
 
 
380
                    raster_matrix[row][col].orig += points->z[0];
 
381
                    raster_matrix[row][col].fo++;
 
382
                }
 
383
 
 
384
                Vect_reset_cats(Cats);
 
385
                Vect_reset_line(points);
 
386
            }
 
387
 
 
388
            for (row = 0; row <= nrows; row++) {
 
389
                for (col = 0; col <= ncols; col++) {
 
390
 
 
391
                    if (raster_matrix[row][col].fc != 0) {
 
392
                        raster_matrix[row][col].clas--;
 
393
                        raster_matrix[row][col].
 
394
                            clas /= raster_matrix[row][col].fc;
 
395
                    }
 
396
 
 
397
                    if (raster_matrix[row][col].fi != 0)
 
398
                        raster_matrix[row][col].
 
399
                            interp /= raster_matrix[row][col].fi;
 
400
 
 
401
                    if (raster_matrix[row][col].fo != 0)
 
402
                        raster_matrix[row][col].
 
403
                            orig /= raster_matrix[row][col].fo;
 
404
                }
 
405
            }
 
406
 
 
407
            /* DOUBLE IMPULSE */
 
408
            Vect_rewind(&First);
 
409
            while (Vect_read_next_line(&First, points_first, Cats_first) > 0) {
 
410
 
 
411
                if ((Vect_point_in_box
 
412
                     (points_first->x[0], points_first->y[0],
 
413
                      points_first->z[0], &elaboration_box)) &&
 
414
                    ((points->x[0] != elaboration_reg.west) ||
 
415
                     (points->x[0] == original_reg.west)) &&
 
416
                    ((points->y[0] != elaboration_reg.north) ||
 
417
                     (points->y[0] == original_reg.north))) {
 
418
 
 
419
                    row =
 
420
                        (int)(Rast_northing_to_row
 
421
                              (points_first->y[0], &elaboration_reg));
 
422
                    col =
 
423
                        (int)(Rast_easting_to_col
 
424
                              (points_first->x[0], &elaboration_reg));
 
425
 
 
426
                    if (fabs
 
427
                        (points_first->z[0] - raster_matrix[row][col].orig) >=
 
428
                        Thres_d)
 
429
                        raster_matrix[row][col].dueImp = DOUBLE_PULSE;
 
430
                }
 
431
                Vect_reset_cats(Cats_first);
 
432
                Vect_reset_line(points_first);
 
433
            }
 
434
 
 
435
            /* REGION GROWING */
 
436
            if (region == TRUE) {
 
437
                G_verbose_message(_("Region Growing"));
 
438
 
 
439
                punti_bordo = G_alloc_matrix(MaxPoints, 3);
 
440
                P = Pvector(0, MaxPoints);
 
441
 
 
442
                colorBordo = 5;
 
443
                ripieno = 6;
 
444
 
 
445
                for (row = 0; row <= nrows; row++) {
 
446
                    G_percent(row, nrows, 2);
 
447
                    for (col = 0; col <= ncols; col++) {
 
448
 
 
449
                        if ((raster_matrix[row][col].clas >= Thres_j) &&
 
450
                            (raster_matrix[row][col].clas < colorBordo)
 
451
                            && (raster_matrix[row][col].fi != 0) &&
 
452
                            (raster_matrix[row][col].dueImp ==
 
453
                             SINGLE_PULSE)) {
 
454
 
 
455
                            /* Selecting a connected Object zone */
 
456
                            ripieno++;
 
457
                            if (ripieno > 10)
 
458
                                ripieno = 6;
 
459
 
 
460
                            /* Selecting points on a connected edge */
 
461
                            for (conta = 0; conta < MaxPoints; conta++) {
 
462
                                punti_bordo[conta][0] = 0;
 
463
                                punti_bordo[conta][1] = 0;
 
464
                                punti_bordo[conta][2] = 0;
 
465
                                P[conta] = punti_bordo[conta];  /* It only makes indexes to be equal, not coord values!! */
 
466
                            }
 
467
 
 
468
                            lungPunti = 0;
 
469
                            lungHull = 0;
 
470
 
 
471
                            regGrow8(elaboration_reg, raster_matrix,
 
472
                                     punti_bordo, &lungPunti, row, col,
 
473
                                     colorBordo, Thres_j, MaxPoints);
 
474
 
 
475
                            /* CONVEX-HULL COMPUTATION */
 
476
                            lungHull = ch2d(P, lungPunti);
 
477
                            cvxHull = G_alloc_matrix(lungHull, 3);
 
478
 
 
479
 
 
480
                            for (xi = 0; xi < lungHull; xi++) {
 
481
                                cvxHull[xi][0] = P[xi][0];
 
482
                                cvxHull[xi][1] = P[xi][1];
 
483
                                cvxHull[xi][2] = P[xi][2];
 
484
                            }
 
485
 
 
486
                            /* Computes the interpoling plane based only on Object points */
 
487
                            altPiano =
 
488
                                pianOriz(punti_bordo, lungPunti, &minNS,
 
489
                                         &minEW, &maxNS, &maxEW,
 
490
                                         raster_matrix, colorBordo);
 
491
 
 
492
                            for (c1 = minNS; c1 <= maxNS; c1++) {
 
493
                                for (c2 = minEW; c2 <= maxEW; c2++) {
 
494
                                    if (checkHull(c1, c2, cvxHull, lungHull)
 
495
                                        == 1) {
 
496
                                        raster_matrix[c1][c2].obj = count_obj;
 
497
 
 
498
                                        if ((raster_matrix[c1][c2].clas ==
 
499
                                             PRE_TERRAIN)
 
500
                                            && (raster_matrix[c1][c2].orig >=
 
501
                                                altPiano) && (lungHull > 3))
 
502
                                            raster_matrix[c1][c2].clas =
 
503
                                                ripieno;
 
504
                                    }
 
505
                                }
 
506
                            }
 
507
                            G_free_matrix(cvxHull);
 
508
                            count_obj++;
 
509
                        }
 
510
                    }
 
511
                }
 
512
                G_free_matrix(punti_bordo);
 
513
                free_Pvector(P, 0, MaxPoints);
 
514
            }
 
515
 
 
516
            /* WRITING THE OUTPUT VECTOR CATEGORIES */
 
517
            Vect_rewind(&In);
 
518
            while (Vect_read_next_line(&In, points, Cats) > 0) {        /* Read every line for buffering points */
 
519
 
 
520
                if ((Vect_point_in_box
 
521
                     (points->x[0], points->y[0], points->z[0],
 
522
                      &elaboration_box)) &&
 
523
                    ((points->x[0] != elaboration_reg.west) ||
 
524
                     (points->x[0] == original_reg.west)) &&
 
525
                    ((points->y[0] != elaboration_reg.north) ||
 
526
                     (points->y[0] == original_reg.north))) {
 
527
 
 
528
                    row =
 
529
                        (int)(Rast_northing_to_row
 
530
                              (points->y[0], &elaboration_reg));
 
531
                    col =
 
532
                        (int)(Rast_easting_to_col
 
533
                              (points->x[0], &elaboration_reg));
 
534
 
 
535
                    if (raster_matrix[row][col].clas == PRE_TERRAIN) {
 
536
                        if (raster_matrix[row][col].dueImp == SINGLE_PULSE)
 
537
                            Vect_cat_set(Cats, F_CLASSIFICATION,
 
538
                                         TERRAIN_SINGLE);
 
539
                        else
 
540
                            Vect_cat_set(Cats, F_CLASSIFICATION,
 
541
                                         TERRAIN_DOUBLE);
 
542
                    }
 
543
                    else {
 
544
                        if (raster_matrix[row][col].dueImp == SINGLE_PULSE)
 
545
                            Vect_cat_set(Cats, F_CLASSIFICATION,
 
546
                                         OBJECT_SINGLE);
 
547
                        else
 
548
                            Vect_cat_set(Cats, F_CLASSIFICATION,
 
549
                                         OBJECT_DOUBLE);
 
550
                    }
 
551
 
 
552
                    Vect_cat_set(Cats, F_COUNTER_OBJ,
 
553
                                 raster_matrix[row][col].obj);
 
554
                    Vect_write_line(&Out, GV_POINT, points, Cats);
 
555
                }
 
556
                Vect_reset_cats(Cats);
 
557
                Vect_reset_line(points);
 
558
            }
 
559
            free_structmatrix(raster_matrix, 0, nrows - 1, 0, ncols - 1);
 
560
        }                       /*! END WHILE; last_column = TRUE */
 
561
    }                           /*! END WHILE; last_row = TRUE */
 
562
 
 
563
    Vect_close(&In);
 
564
    Vect_close(&First);
 
565
    Vect_close(&Out);
 
566
 
 
567
    db_close_database_shutdown_driver(driver);
 
568
 
 
569
    G_done_msg(" ");
 
570
    exit(EXIT_SUCCESS);
 
571
}