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_LINE *Line;
00336 struct Plus_head *plus;
00337
00338 G_debug(3, "Vect_attach_centroids ()");
00339
00340 plus = &(Map->plus);
00341
00342 if (first) {
00343 List = Vect_new_list();
00344 first = 0;
00345 }
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
00371
00372
00373 Vect_select_lines_by_box(Map, box, GV_CENTROID, List);
00374 G_debug(3, " number of centroids to reattach = %d", List->n_values);
00375 for (i = 0; i < List->n_values; i++) {
00376 int orig_area;
00377
00378 centr = List->value[i];
00379 Line = plus->Line[centr];
00380
00381
00382
00383
00384
00385
00386 if (Line->left > 0)
00387 continue;
00388
00389 orig_area = Line->left;
00390
00391 sel_area = Vect_find_area(Map, Line->E, Line->N);
00392 G_debug(3, " centroid %d is in area %d", centr, sel_area);
00393 if (sel_area > 0) {
00394 Area = plus->Area[sel_area];
00395 if (Area->centroid == 0) {
00396 G_debug(3, " first centroid -> attach to area");
00397 Area->centroid = centr;
00398 Line->left = sel_area;
00399
00400 if (sel_area != orig_area && plus->do_uplist)
00401 dig_line_add_updated(plus, centr);
00402 }
00403 else if (Area->centroid != centr) {
00404
00405
00406 G_debug(3, " duplicate centroid -> do not attach to area");
00407 Line->left = -sel_area;
00408
00409 if (-sel_area != orig_area && plus->do_uplist)
00410 dig_line_add_updated(plus, centr);
00411 }
00412 }
00413
00414 if (sel_area != orig_area && plus->do_uplist)
00415 dig_line_add_updated(plus, centr);
00416 }
00417
00418 return 0;
00419 }
00420
00430 int Vect_build_nat(struct Map_info *Map, int build)
00431 {
00432 struct Plus_head *plus;
00433 int i, s, type, lineid;
00434 long offset;
00435 int side, line, area;
00436 struct line_pnts *Points, *APoints;
00437 struct line_cats *Cats;
00438 P_LINE *Line;
00439 P_NODE *Node;
00440 P_AREA *Area;
00441 BOUND_BOX box;
00442 struct ilist *List;
00443
00444 G_debug(3, "Vect_build_nat() build = %d", build);
00445
00446 plus = &(Map->plus);
00447
00448 if (build == plus->built)
00449 return 1;
00450
00451
00452 if (build < plus->built) {
00453
00454
00455 if (plus->built >= GV_BUILD_CENTROIDS && build < GV_BUILD_CENTROIDS) {
00456
00457 int nlines = Vect_get_num_lines(Map);
00458
00459 for (line = 1; line <= nlines; line++) {
00460 Line = plus->Line[line];
00461 if (Line && Line->type == GV_CENTROID)
00462 Line->left = 0;
00463 }
00464 dig_free_plus_areas(plus);
00465 dig_spidx_free_areas(plus);
00466 dig_free_plus_isles(plus);
00467 dig_spidx_free_isles(plus);
00468 }
00469
00470
00471 if (plus->built >= GV_BUILD_AREAS && build < GV_BUILD_AREAS) {
00472
00473 int nlines = Vect_get_num_lines(Map);
00474
00475 for (line = 1; line <= nlines; line++) {
00476 Line = plus->Line[line];
00477 if (Line && Line->type == GV_BOUNDARY) {
00478 Line->left = 0;
00479 Line->right = 0;
00480 }
00481 }
00482 dig_free_plus_areas(plus);
00483 dig_spidx_free_areas(plus);
00484 dig_free_plus_isles(plus);
00485 dig_spidx_free_isles(plus);
00486 }
00487 if (plus->built >= GV_BUILD_BASE && build < GV_BUILD_BASE) {
00488 dig_free_plus_nodes(plus);
00489 dig_spidx_free_nodes(plus);
00490 dig_free_plus_lines(plus);
00491 dig_spidx_free_lines(plus);
00492 }
00493
00494 plus->built = build;
00495 return 1;
00496 }
00497
00498 Points = Vect_new_line_struct();
00499 APoints = Vect_new_line_struct();
00500 Cats = Vect_new_cats_struct();
00501 List = Vect_new_list();
00502
00503 if (plus->built < GV_BUILD_BASE) {
00504 int npoints, format;
00505
00506 format = G_info_format();
00507
00508
00509
00510
00511
00512
00513
00514
00515 Vect_rewind(Map);
00516 G_message(_("Registering primitives..."));
00517 i = 1;
00518 npoints = 0;
00519 while (1) {
00520
00521 type = Vect_read_next_line(Map, Points, Cats);
00522
00523
00524 if (type == -1) {
00525 G_warning(_("Unable to read vector map"));
00526 return 0;
00527 }
00528 else if (type == -2) {
00529 break;
00530 }
00531
00532 npoints += Points->n_points;
00533
00534 offset = Map->head.last_offset;
00535
00536 G_debug(3, "Register line: offset = %ld", offset);
00537 lineid = dig_add_line(plus, type, Points, offset);
00538 dig_line_box(Points, &box);
00539 if (lineid == 1)
00540 Vect_box_copy(&(plus->box), &box);
00541 else
00542 Vect_box_extend(&(plus->box), &box);
00543
00544
00545 if (build == GV_BUILD_ALL) {
00546 int c;
00547
00548 for (c = 0; c < Cats->n_cats; c++) {
00549 dig_cidx_add_cat(plus, Cats->field[c], Cats->cat[c],
00550 lineid, type);
00551 }
00552 if (Cats->n_cats == 0)
00553 dig_cidx_add_cat(plus, 0, 0, lineid, type);
00554 }
00555
00556 if (G_verbose() > G_verbose_min() && i % 1000 == 0) {
00557 if (format == G_INFO_FORMAT_PLAIN)
00558 fprintf(stderr, "%d..", i);
00559 else
00560 fprintf(stderr, "%7d\b\b\b\b\b\b\b", i);
00561 }
00562
00563 i++;
00564 }
00565
00566 if ( (G_verbose() > G_verbose_min() ) && format != G_INFO_FORMAT_PLAIN )
00567 fprintf(stderr, "\r");
00568
00569 G_message(_("%d primitives registered"), plus->n_lines);
00570 G_message(_("%d vertices registered"), npoints);
00571
00572 plus->built = GV_BUILD_BASE;
00573 }
00574
00575 if (build < GV_BUILD_AREAS)
00576 return 1;
00577
00578 if (plus->built < GV_BUILD_AREAS) {
00579
00580
00581 G_message(_("Building areas..."));
00582 for (i = 1; i <= plus->n_lines; i++) {
00583 G_percent(i, plus->n_lines, 1);
00584
00585
00586 if (plus->Line[i] == NULL) {
00587 continue;
00588 }
00589 Line = plus->Line[i];
00590 if (Line->type != GV_BOUNDARY) {
00591 continue;
00592 }
00593
00594 for (s = 0; s < 2; s++) {
00595 if (s == 0)
00596 side = GV_LEFT;
00597 else
00598 side = GV_RIGHT;
00599
00600 G_debug(3, "Build area for line = %d, side = %d", i, side);
00601 Vect_build_line_area(Map, i, side);
00602 }
00603 }
00604 G_message(_("%d areas built"), plus->n_areas);
00605 G_message(_("%d isles built"), plus->n_isles);
00606 plus->built = GV_BUILD_AREAS;
00607 }
00608
00609 if (build < GV_BUILD_ATTACH_ISLES)
00610 return 1;
00611
00612
00613 if (plus->built < GV_BUILD_ATTACH_ISLES) {
00614 G_message(_("Attaching islands..."));
00615 for (i = 1; i <= plus->n_isles; i++) {
00616 G_percent(i, plus->n_isles, 1);
00617 Vect_attach_isle(Map, i);
00618 }
00619 plus->built = GV_BUILD_ATTACH_ISLES;
00620 }
00621
00622 if (build < GV_BUILD_CENTROIDS)
00623 return 1;
00624
00625
00626 if (plus->built < GV_BUILD_CENTROIDS) {
00627 int nlines;
00628
00629 G_message(_("Attaching centroids..."));
00630
00631 nlines = Vect_get_num_lines(Map);
00632 for (line = 1; line <= nlines; line++) {
00633 G_percent(line, nlines, 1);
00634
00635 Line = plus->Line[line];
00636 if (!Line)
00637 continue;
00638
00639 if (Line->type != GV_CENTROID)
00640 continue;
00641
00642 Node = plus->Node[Line->N1];
00643
00644 area = Vect_find_area(Map, Node->x, Node->y);
00645
00646 if (area > 0) {
00647 G_debug(3, "Centroid (line=%d) in area %d", line, area);
00648
00649 Area = plus->Area[area];
00650
00651 if (Area->centroid == 0) {
00652 Area->centroid = line;
00653 Line->left = area;
00654 }
00655 else {
00656 Line->left = -area;
00657 }
00658 }
00659 }
00660 plus->built = GV_BUILD_CENTROIDS;
00661 }
00662
00663
00664 for (area = 1; area <= plus->n_areas; area++) {
00665 int c;
00666
00667 if (plus->Area[area] == NULL)
00668 continue;
00669
00670 if (plus->Area[area]->centroid > 0) {
00671 Vect_read_line(Map, NULL, Cats, plus->Area[area]->centroid);
00672
00673 for (c = 0; c < Cats->n_cats; c++) {
00674 dig_cidx_add_cat(plus, Cats->field[c], Cats->cat[c], area,
00675 GV_AREA);
00676 }
00677 }
00678
00679 if (plus->Area[area]->centroid == 0 || Cats->n_cats == 0)
00680 dig_cidx_add_cat(plus, 0, 0, area, GV_AREA);
00681 }
00682
00683 return 1;
00684 }