00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <string.h>
00019 #include <stdlib.h>
00020 #include <stdio.h>
00021 #include <grass/glocale.h>
00022 #include <grass/gis.h>
00023 #include <grass/Vect.h>
00024
00025 extern FILE *Msgout;
00026 extern int prnmsg ( char *msg, ...) ;
00027
00034 int
00035 Vect_build_line_area ( struct Map_info *Map, int iline, int side )
00036 {
00037 int j, area, isle, n_lines, line, type, direction;
00038 static int first = 1;
00039 long offset;
00040 struct Plus_head *plus ;
00041 P_LINE *BLine;
00042 static struct line_pnts *Points, *APoints;
00043 plus_t *lines;
00044 double area_size;
00045
00046 plus = &(Map->plus);
00047
00048 G_debug ( 3, "Vect_build_line_area() line = %d, side = %d", iline, side );
00049
00050 if ( first ) {
00051 Points = Vect_new_line_struct ();
00052 APoints = Vect_new_line_struct ();
00053 first = 0;
00054 }
00055
00056 area = dig_line_get_area (plus, iline, side);
00057 if ( area != 0 ) {
00058 G_debug ( 3, " area/isle = %d -> skip", area );
00059 return 0;
00060 }
00061
00062 n_lines = dig_build_area_with_line (plus, iline, side, &lines);
00063 G_debug ( 3, " n_lines = %d", n_lines );
00064 if ( n_lines < 1 ) { return 0; }
00065
00066
00067 Vect_reset_line ( APoints );
00068 for (j = 0; j < n_lines; j++){
00069 line = abs(lines[j]);
00070 BLine = plus->Line[line];
00071 offset = BLine->offset;
00072 G_debug ( 3, " line[%d] = %d, offset = %d", j, line, offset );
00073 type = Vect_read_line (Map, Points, NULL, line );
00074 if ( lines[j] > 0 ) direction = GV_FORWARD;
00075 else direction = GV_BACKWARD;
00076 Vect_append_points ( APoints, Points, direction);
00077 }
00078
00079 dig_find_area_poly (APoints, &area_size);
00080 G_debug ( 3, " area/isle size = %f", area_size );
00081
00082 if (area_size > 0) {
00083
00084 area = dig_add_area (plus, n_lines, lines);
00085 if ( area == -1 ) {
00086 Vect_close ( Map );
00087 G_fatal_error ( _("Cannot add area (map closed, topo saved)") );
00088 }
00089 G_debug ( 3, " -> area %d", area );
00090 return area;
00091 } else if (area_size < 0) {
00092 isle = dig_add_isle (plus, n_lines, lines);
00093 if ( isle == -1 ) {
00094 Vect_close ( Map );
00095 G_fatal_error ( _("Cannot add isle (map closed, topo saved)") );
00096 }
00097 G_debug ( 3, " -> isle %d", isle );
00098 return -isle;
00099 } else {
00100
00101
00102
00103
00104 G_warning (_("Area of size = 0.0 ignored"));
00105 }
00106 return 0;
00107 }
00108
00115 int
00116 Vect_isle_find_area ( struct Map_info *Map, int isle )
00117 {
00118 int j, line, node, sel_area, first, area, poly;
00119 static int first_call = 1;
00120 struct Plus_head *plus ;
00121 P_LINE *Line;
00122 P_NODE *Node;
00123 P_ISLE *Isle;
00124 P_AREA *Area;
00125 double size, cur_size;
00126 BOUND_BOX box, abox;
00127 static struct ilist *List;
00128 static struct line_pnts *APoints;
00129
00130
00131
00132
00133
00134 G_debug ( 3, "Vect_isle_find_area () island = %d", isle );
00135 plus = &(Map->plus);
00136
00137 if ( plus->Isle[isle] == NULL ) {
00138 G_warning ("Request to find area outside nonexisting isle");
00139 return 0;
00140 }
00141
00142 if ( first_call ) {
00143 List = Vect_new_list ();
00144 APoints = Vect_new_line_struct ();
00145 first_call = 0;
00146 }
00147
00148 Isle = plus->Isle[isle];
00149 line = abs(Isle->lines[0]);
00150 Line = plus->Line[line];
00151 node = Line->N1;
00152 Node = plus->Node[node];
00153
00154
00155 box.E = Node->x; box.W = Node->x; box.N = Node->y; box.S = Node->y;
00156 box.T = PORT_DOUBLE_MAX; box.B = -PORT_DOUBLE_MAX;
00157 Vect_select_areas_by_box (Map, &box, List);
00158 G_debug ( 3, "%d areas overlap island boundary point", List->n_values );
00159
00160 sel_area = 0; cur_size = -1;
00161 first = 1;
00162 Vect_get_isle_box ( Map, isle, &box);
00163 for (j = 0; j < List->n_values; j++) {
00164 area = List->value[j];
00165 G_debug ( 3, "area = %d", area );
00166
00167 Area = plus->Area[area];
00168
00169
00170 if ( abs(Isle->lines[0]) == abs(Area->lines[0]) ) {
00171 G_debug ( 3, " area inside isolated isle" );
00172 continue;
00173 }
00174
00175
00176
00177
00178
00179
00180
00181 Vect_get_area_box ( Map, area, &abox);
00182 if ( box.E > abox.E || box.W < abox.W || box.N > abox.N || box.S < abox.S ) {
00183 G_debug ( 3, " isle not completely inside area box" );
00184 continue;
00185 }
00186
00187 poly = Vect_point_in_area_outer_ring ( Node->x, Node->y, Map, area);
00188 G_debug ( 3, " poly = %d", poly );
00189
00190 if ( poly == 1) {
00191
00192
00193 if ( sel_area == 0 ) {
00194 sel_area = area;
00195 } else {
00196 if ( cur_size < 0 ) {
00197
00198 Vect_get_area_points (Map, sel_area, APoints);
00199 G_begin_polygon_area_calculations();
00200 cur_size = G_area_of_polygon(APoints->x, APoints->y, APoints->n_points);
00201 G_debug ( 3, " first area size = %f (n points = %d)", cur_size, APoints->n_points );
00202
00203 }
00204
00205 Vect_get_area_points (Map, area, APoints);
00206 size = G_area_of_polygon(APoints->x, APoints->y, APoints->n_points);
00207 G_debug ( 3, " area size = %f (n points = %d)", cur_size, APoints->n_points );
00208
00209 if ( size < cur_size ) {
00210 sel_area = area;
00211 cur_size = size;
00212 }
00213 }
00214 G_debug ( 3, "sel_area = %d cur_size = %f", sel_area, cur_size );
00215 }
00216 }
00217 if ( sel_area > 0 ) {
00218 G_debug ( 3, "Island %d in area %d", isle, sel_area );
00219 } else {
00220 G_debug ( 3, "Island %d is not in area", isle );
00221 }
00222
00223 return sel_area;
00224 }
00225
00232 int
00233 Vect_attach_isle ( struct Map_info *Map, int isle )
00234 {
00235 int sel_area;
00236 P_ISLE *Isle;
00237 struct Plus_head *plus ;
00238
00239
00240
00241 G_debug ( 3, "Vect_attach_isle (): isle = %d", isle);
00242
00243 plus = &(Map->plus);
00244
00245 sel_area = Vect_isle_find_area ( Map, isle );
00246 G_debug ( 3, " isle = %d -> area outside = %d", isle, sel_area);
00247 if ( sel_area > 0 ) {
00248 Isle = plus->Isle[isle];
00249 if ( Isle->area > 0 ) {
00250 G_debug (3, "Attempt to attach isle %d to more areas (=>topology is not clean)");
00251 } else {
00252 Isle->area = sel_area;
00253 dig_area_add_isle ( plus, sel_area, isle );
00254 }
00255 }
00256 return 0;
00257 }
00258
00265 int
00266 Vect_attach_isles ( struct Map_info *Map, BOUND_BOX *box )
00267 {
00268 int i, isle;
00269 static int first = 1;
00270 static struct ilist *List;
00271 struct Plus_head *plus ;
00272
00273 G_debug ( 3, "Vect_attach_isles ()");
00274
00275 plus = &(Map->plus);
00276
00277 if ( first ) {
00278 List = Vect_new_list ();
00279 first = 0;
00280 }
00281
00282 Vect_select_isles_by_box ( Map, box, List);
00283 G_debug ( 3, " number of isles to attach = %d", List->n_values);
00284
00285 for (i = 0; i < List->n_values; i++) {
00286 isle = List->value[i];
00287 Vect_attach_isle ( Map, isle );
00288 }
00289 return 0;
00290 }
00291
00298 int
00299 Vect_attach_centroids ( struct Map_info *Map, BOUND_BOX *box )
00300 {
00301 int i, sel_area, centr;
00302 static int first = 1;
00303 static struct ilist *List;
00304 P_AREA *Area;
00305 P_NODE *Node;
00306 P_LINE *Line;
00307 struct Plus_head *plus ;
00308
00309 G_debug ( 3, "Vect_attach_centroids ()");
00310
00311 plus = &(Map->plus);
00312
00313 if ( first ) {
00314 List = Vect_new_list ();
00315 first = 0;
00316 }
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340 Vect_select_lines_by_box ( Map, box, GV_CENTROID, List);
00341 G_debug ( 3, " number of centroids to reattach = %d", List->n_values);
00342 for (i = 0; i < List->n_values; i++) {
00343 int orig_area;
00344 centr = List->value[i];
00345 Line = plus->Line[centr];
00346 Node = plus->Node[Line->N1];
00347
00348
00349 orig_area = Line->left;
00350 if ( orig_area > 0 ) {
00351 if ( plus->Area[orig_area] != NULL ) {
00352 plus->Area[orig_area]->centroid = 0;
00353 }
00354 }
00355 Line->left = 0;
00356
00357 sel_area = Vect_find_area ( Map, Node->x, Node->y );
00358 G_debug ( 3, " centroid %d is in area %d", centr, sel_area);
00359 if ( sel_area > 0 ) {
00360 Area = plus->Area[sel_area];
00361 if ( Area->centroid == 0 ) {
00362 G_debug ( 3, " first centroid -> attach to area");
00363 Area->centroid = centr;
00364 Line->left = sel_area;
00365 } else if ( Area->centroid != centr ) {
00366
00367
00368 G_debug ( 3, " duplicate centroid -> do not attach to area");
00369 Line->left = -sel_area;
00370 }
00371 }
00372
00373 if ( sel_area != orig_area && plus->do_uplist ) dig_line_add_updated ( plus, centr );
00374 }
00375
00376 return 0;
00377 }
00378
00385 int
00386 Vect_build_nat ( struct Map_info *Map, int build, FILE *msgout )
00387 {
00388 struct Plus_head *plus ;
00389 int i, j, s, type, lineid;
00390 long offset;
00391 int side, line, area;
00392 struct line_pnts *Points, *APoints;
00393 struct line_cats *Cats;
00394 P_LINE *Line;
00395 P_NODE *Node;
00396 P_AREA *Area;
00397 BOUND_BOX box;
00398 struct ilist *List;
00399
00400 G_debug (3, "Vect_build_nat() build = %d", build);
00401
00402 plus = &(Map->plus);
00403 Msgout = msgout;
00404
00405 if ( build == plus->built ) return 1;
00406
00407
00408 if ( build < plus->built ) {
00409
00410
00411 if ( plus->built >= GV_BUILD_CENTROIDS && build < GV_BUILD_CENTROIDS) {
00412
00413 int nlines = Vect_get_num_lines (Map);
00414 for ( line = 1; line <= nlines; line++ ) {
00415 Line = plus->Line[line];
00416 if ( Line && Line->type == GV_CENTROID )
00417 Line->left = 0;
00418 }
00419 dig_free_plus_areas (plus);
00420 dig_spidx_free_areas (plus);
00421 dig_free_plus_isles (plus);
00422 dig_spidx_free_isles (plus);
00423 }
00424
00425
00426 if ( plus->built >= GV_BUILD_AREAS && build < GV_BUILD_AREAS) {
00427
00428 int nlines = Vect_get_num_lines (Map);
00429 for ( line = 1; line <= nlines; line++ ) {
00430 Line = plus->Line[line];
00431 if ( Line && Line->type == GV_BOUNDARY ) {
00432 Line->left = 0;
00433 Line->right = 0;
00434 }
00435 }
00436 dig_free_plus_areas (plus);
00437 dig_spidx_free_areas (plus);
00438 dig_free_plus_isles (plus);
00439 dig_spidx_free_isles (plus);
00440 }
00441 if ( plus->built >= GV_BUILD_BASE && build < GV_BUILD_BASE) {
00442 dig_free_plus_nodes (plus);
00443 dig_spidx_free_nodes (plus);
00444 dig_free_plus_lines (plus);
00445 dig_spidx_free_lines (plus);
00446 }
00447
00448 plus->built = build;
00449 return 1;
00450 }
00451
00452 Points = Vect_new_line_struct ();
00453 APoints = Vect_new_line_struct ();
00454 Cats = Vect_new_cats_struct ();
00455 List = Vect_new_list ();
00456
00457 if ( plus->built < GV_BUILD_BASE ) {
00458
00459
00460
00461
00462
00463
00464
00465 Vect_rewind ( Map );
00466 prnmsg (_("Registering lines: "));
00467 i = 1; j = 1;
00468 while ( 1 ) {
00469
00470 type = Vect_read_next_line (Map, Points, Cats);
00471
00472 if ( type == -1 ) {
00473 fprintf (stderr, "\nERROR: vector file - can't read\n" );
00474 return 0;
00475 } else if ( type == -2 ) {
00476 break;
00477 }
00478
00479 offset = Map->head.last_offset;
00480
00481 G_debug ( 3, "Register line: offset = %ld", offset );
00482 lineid = dig_add_line ( plus, type, Points, offset );
00483 dig_line_box ( Points, &box );
00484 if ( lineid == 1 )
00485 Vect_box_copy (&(plus->box), &box);
00486 else
00487 Vect_box_extend (&(plus->box), &box);
00488
00489
00490 if ( build == GV_BUILD_ALL ) {
00491 int c;
00492
00493 for ( c = 0; c < Cats->n_cats; c++ ) {
00494 dig_cidx_add_cat (plus, Cats->field[c], Cats->cat[c], lineid, type);
00495 }
00496 if ( Cats->n_cats == 0 )
00497 dig_cidx_add_cat (plus, 0, 0, lineid, type);
00498 }
00499
00500
00501 if ( i == 1000 ) {
00502 prnmsg ("%7d\b\b\b\b\b\b\b", j);
00503 i = 0;
00504 }
00505 i++; j++;
00506 }
00507 prnmsg (_("\r%d primitives registered \n"), plus->n_lines);
00508
00509 plus->built = GV_BUILD_BASE;
00510 }
00511
00512 if ( build < GV_BUILD_AREAS ) return 1;
00513
00514 if ( plus->built < GV_BUILD_AREAS ) {
00515
00516
00517 prnmsg ("Building areas: ");
00518 for (i = 1; i <= plus->n_lines; i++) {
00519 G_percent2 ( i, plus->n_lines, 1, msgout );
00520
00521
00522 if ( plus->Line[i] == NULL ) { continue; }
00523 Line = plus->Line[i];
00524 if ( Line->type != GV_BOUNDARY ) { continue; }
00525
00526 for (s = 0; s < 2; s++) {
00527 if ( s == 0 ) side = GV_LEFT;
00528 else side = GV_RIGHT;
00529
00530 G_debug ( 3, "Build area for line = %d, side = %d", i, side );
00531 Vect_build_line_area ( Map, i, side );
00532 }
00533 }
00534 prnmsg (_("\r%d areas built \n%d isles built\n"), plus->n_areas, plus->n_isles );
00535 plus->built = GV_BUILD_AREAS;
00536 }
00537
00538 if ( build < GV_BUILD_ATTACH_ISLES ) return 1;
00539
00540
00541 if ( plus->built < GV_BUILD_ATTACH_ISLES ) {
00542 prnmsg (_("Attaching islands: "));
00543 for (i = 1; i <= plus->n_isles; i++) {
00544 G_percent2 ( i, plus->n_isles, 1, msgout );
00545 Vect_attach_isle ( Map, i ) ;
00546 }
00547 if (i==1) prnmsg ("\n");
00548 plus->built = GV_BUILD_ATTACH_ISLES;
00549 }
00550
00551 if ( build < GV_BUILD_CENTROIDS ) return 1;
00552
00553
00554 if ( plus->built < GV_BUILD_CENTROIDS ) {
00555 int nlines;
00556
00557 prnmsg (_("Attaching centroids: "));
00558
00559 nlines = Vect_get_num_lines (Map);
00560 for ( line = 1; line <= nlines; line++ ) {
00561 G_percent2 ( line, nlines, 1, msgout );
00562
00563 Line = plus->Line[line];
00564 if ( !Line ) continue;
00565
00566 if ( Line->type != GV_CENTROID ) continue;
00567
00568 Node = plus->Node[Line->N1];
00569
00570 area = Vect_find_area ( Map, Node->x, Node->y );
00571
00572 if ( area > 0 ) {
00573 G_debug ( 3, "Centroid (line=%d) in area %d", line, area );
00574
00575 Area = plus->Area[area];
00576
00577 if ( Area->centroid == 0 ) {
00578 Area->centroid = line;
00579 Line->left = area;
00580 } else {
00581 Line->left = -area;
00582 }
00583 }
00584 }
00585 plus->built = GV_BUILD_CENTROIDS;
00586 }
00587
00588
00589 for ( area = 1; area <= plus->n_areas; area++ ) {
00590 int c;
00591
00592 if ( plus->Area[area] == NULL ) continue;
00593
00594 if ( plus->Area[area]->centroid > 0 ) {
00595 Vect_read_line (Map, NULL, Cats, plus->Area[area]->centroid );
00596
00597 for ( c = 0; c < Cats->n_cats; c++ ) {
00598 dig_cidx_add_cat (plus, Cats->field[c], Cats->cat[c], area, GV_AREA);
00599 }
00600 }
00601
00602 if ( plus->Area[area]->centroid == 0 || Cats->n_cats == 0 )
00603 dig_cidx_add_cat (plus, 0, 0, area, GV_AREA);
00604 }
00605
00606 return 1;
00607 }
00608
00609