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