2
/****************************************************************
4
* MODULE: v.lidar.growing *
6
* AUTHOR(S): Roberto Antolin & Gonzalo Moreno *
8
* PURPOSE: Building contour determination and Region *
9
* Growing algorithm for determining the building *
12
* COPYRIGHT: (C) 2006 by Politecnico di Milano - *
13
* Polo Regionale di Como *
15
* This program is free software under the *
16
* GNU General Public License (>=v2). *
17
* Read the file COPYING that comes with GRASS *
20
****************************************************************/
28
/* GLOBAL DEFINITIONS */
29
int nsply, nsplx, count_obj;
32
/*--------------------------------------------------------------------------------*/
33
int main(int argc, char *argv[])
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;
43
int clas, region = TRUE;
45
double Thres_j, Thres_d, ew_resol, ns_resol;
46
double minNS, minEW, maxNS, maxEW;
48
char buf[1024], table_name[GNAME_MAX];
49
char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
51
int colorBordo, ripieno, conta, lungPunti, lungHull, xi, c1, c2;
53
extern double **P, **cvxHull, **punti_bordo;
55
/* Struct declarations */
56
struct Cell_head elaboration_reg, original_reg;
57
struct element_grow **raster_matrix;
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;
63
struct line_pnts *points, *points_first;
64
struct line_cats *Cats, *Cats_first;
66
struct field_info *field;
72
/*------------------------------------------------------------------------------------------*/
73
/* Options' declaration */ ;
74
module = G_define_module();
75
G_add_keyword(_("vector"));
76
G_add_keyword(_("LIDAR"));
78
_("Building contour determination and Region Growing "
79
"algorithm for determining the building inside");
81
in_opt = G_define_standard_option(G_OPT_V_INPUT);
83
_("Input vector (v.lidar.edgedetection output)");
85
out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
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");
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";
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";
113
if (G_parser(argc, argv))
116
Thres_j = atof(Thres_j_opt->answer);
117
Thres_d = atof(Thres_d_opt->answer);
121
/* Open input vector */
122
Vect_check_input_output_name(in_opt->answer, out_opt->answer,
125
if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL) {
126
G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
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);
134
sprintf(table_name, "%s_edge_Interpolation", in_opt->answer);
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);
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);
144
/* Open output vector */
145
if (0 > Vect_open_new(&Out, out_opt->answer, WITH_Z)) {
151
/* Copy vector Head File */
152
Vect_copy_head_data(&In, &Out);
153
Vect_hist_copy(&In, &Out);
154
Vect_hist_command(&Out);
156
/* Starting driver and open db for edgedetection interpolation table */
157
field = Vect_get_field(&In, F_INTERPOLATION);
159
G_fatal_error (_("Cannot read field info")); */
161
driver = db_start_driver_open_database(field->driver, field->database);
163
G_fatal_error(_("No database connection for driver <%s> is defined. Run db.connect."),
166
/* is this the right place to open the cursor ??? */
168
db_init_string(&sql);
169
db_zero_string(&sql);
171
sprintf(buf, "SELECT Interp,ID FROM %s", table_name);
172
G_debug(1, "buf: %s", buf);
173
db_append_string(&sql, buf);
175
if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
176
G_fatal_error(_("Unable to open table <%s>"), table_name);
180
/* no topology, get number of lines in input vector */
182
points = Vect_new_line_struct();
183
Cats = Vect_new_cats_struct();
185
while (Vect_read_next_line(&In, points, Cats) > 0) {
190
/* no topology, get number of lines in first pulse input vector */
192
points_first = Vect_new_line_struct();
193
Cats_first = Vect_new_cats_struct();
195
while (Vect_read_next_line(&First, points_first, Cats_first) > 0) {
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);
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;
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;
214
if (nsubregion_col < 0)
216
if (nsubregion_row < 0)
219
nsubregions = nsubregion_row * nsubregion_col;
221
/* Subdividing and working with tiles */
222
elaboration_reg.south = original_reg.north;
225
while (last_row == FALSE) { /* For each strip of LATO rows */
227
elaboration_reg.north = elaboration_reg.south;
229
if (elaboration_reg.north > original_reg.north) /* First row */
230
elaboration_reg.north = original_reg.north;
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;
238
elaboration_reg.east = original_reg.west;
241
while (last_column == FALSE) { /* For each strip of LATO columns */
242
struct bound_box elaboration_box;
246
G_message(_("subregion %d of %d"), subregion, nsubregions);
248
elaboration_reg.west = elaboration_reg.east;
249
if (elaboration_reg.west < original_reg.west) /* First column */
250
elaboration_reg.west = original_reg.west;
252
elaboration_reg.east = elaboration_reg.west + LATO * ew_resol;
254
if (elaboration_reg.east >= original_reg.east) { /* Last column */
255
elaboration_reg.east = original_reg.east;
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;
267
G_debug(1, _("Rows = %d"), nrows);
268
G_debug(1, _("Columns = %d"), ncols);
270
raster_matrix = structMatrix(0, nrows, 0, ncols);
271
MaxPoints = nrows * ncols;
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;
288
G_verbose_message(_("read points in input vector"));
289
Vect_region_box(&elaboration_reg, &elaboration_box);
292
while (Vect_read_next_line(&In, points, Cats) > 0) {
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))) {
304
(int)(Rast_northing_to_row
305
(points->y[0], &elaboration_reg));
307
(int)(Rast_easting_to_col
308
(points->x[0], &elaboration_reg));
311
/* TODO: make sure the current db_fetch() usage works */
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);
318
if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK)
319
G_fatal_error(_("Unable to open table <%s>"), table_name);
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);
326
Z_Interp_col = db_get_table_column(table, 1);
328
if (db_sqltype_to_Ctype(db_get_column_sqltype(Z_Interp_col)) ==
330
Z_Interp_value = db_get_column_value(Z_Interp_col);
334
Z_interp = db_get_value_double(Z_Interp_value);
337
db_close_cursor(&cursor);
338
db_free_string(&sql);
342
if (db_fetch(&cursor, DB_NEXT, &more) != DB_OK ||
345
dbColumn *Z_Interp_col, *ID_col;
346
dbValue *Z_Interp_value, *ID_value;
348
table = db_get_cursor_table(&cursor);
350
ID_col = db_get_table_column(table, 1);
351
if (db_sqltype_to_Ctype(db_get_column_sqltype(ID_col))
353
ID_value = db_get_column_value(ID_col);
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)) ==
363
db_get_column_value(Z_Interp_col);
366
Z_interp = db_get_value_double(Z_Interp_value);
371
raster_matrix[row][col].interp += Z_interp;
372
raster_matrix[row][col].fi++;
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++;
380
raster_matrix[row][col].orig += points->z[0];
381
raster_matrix[row][col].fo++;
384
Vect_reset_cats(Cats);
385
Vect_reset_line(points);
388
for (row = 0; row <= nrows; row++) {
389
for (col = 0; col <= ncols; col++) {
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;
397
if (raster_matrix[row][col].fi != 0)
398
raster_matrix[row][col].
399
interp /= raster_matrix[row][col].fi;
401
if (raster_matrix[row][col].fo != 0)
402
raster_matrix[row][col].
403
orig /= raster_matrix[row][col].fo;
409
while (Vect_read_next_line(&First, points_first, Cats_first) > 0) {
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))) {
420
(int)(Rast_northing_to_row
421
(points_first->y[0], &elaboration_reg));
423
(int)(Rast_easting_to_col
424
(points_first->x[0], &elaboration_reg));
427
(points_first->z[0] - raster_matrix[row][col].orig) >=
429
raster_matrix[row][col].dueImp = DOUBLE_PULSE;
431
Vect_reset_cats(Cats_first);
432
Vect_reset_line(points_first);
436
if (region == TRUE) {
437
G_verbose_message(_("Region Growing"));
439
punti_bordo = G_alloc_matrix(MaxPoints, 3);
440
P = Pvector(0, MaxPoints);
445
for (row = 0; row <= nrows; row++) {
446
G_percent(row, nrows, 2);
447
for (col = 0; col <= ncols; col++) {
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 ==
455
/* Selecting a connected Object zone */
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!! */
471
regGrow8(elaboration_reg, raster_matrix,
472
punti_bordo, &lungPunti, row, col,
473
colorBordo, Thres_j, MaxPoints);
475
/* CONVEX-HULL COMPUTATION */
476
lungHull = ch2d(P, lungPunti);
477
cvxHull = G_alloc_matrix(lungHull, 3);
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];
486
/* Computes the interpoling plane based only on Object points */
488
pianOriz(punti_bordo, lungPunti, &minNS,
489
&minEW, &maxNS, &maxEW,
490
raster_matrix, colorBordo);
492
for (c1 = minNS; c1 <= maxNS; c1++) {
493
for (c2 = minEW; c2 <= maxEW; c2++) {
494
if (checkHull(c1, c2, cvxHull, lungHull)
496
raster_matrix[c1][c2].obj = count_obj;
498
if ((raster_matrix[c1][c2].clas ==
500
&& (raster_matrix[c1][c2].orig >=
501
altPiano) && (lungHull > 3))
502
raster_matrix[c1][c2].clas =
507
G_free_matrix(cvxHull);
512
G_free_matrix(punti_bordo);
513
free_Pvector(P, 0, MaxPoints);
516
/* WRITING THE OUTPUT VECTOR CATEGORIES */
518
while (Vect_read_next_line(&In, points, Cats) > 0) { /* Read every line for buffering points */
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))) {
529
(int)(Rast_northing_to_row
530
(points->y[0], &elaboration_reg));
532
(int)(Rast_easting_to_col
533
(points->x[0], &elaboration_reg));
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,
540
Vect_cat_set(Cats, F_CLASSIFICATION,
544
if (raster_matrix[row][col].dueImp == SINGLE_PULSE)
545
Vect_cat_set(Cats, F_CLASSIFICATION,
548
Vect_cat_set(Cats, F_CLASSIFICATION,
552
Vect_cat_set(Cats, F_COUNTER_OBJ,
553
raster_matrix[row][col].obj);
554
Vect_write_line(&Out, GV_POINT, points, Cats);
556
Vect_reset_cats(Cats);
557
Vect_reset_line(points);
559
free_structmatrix(raster_matrix, 0, nrows - 1, 0, ncols - 1);
560
} /*! END WHILE; last_column = TRUE */
561
} /*! END WHILE; last_row = TRUE */
567
db_close_database_shutdown_driver(driver);