[GRASS-SVN] r49179 - grass/trunk/lib/vector/Vlib
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Nov 11 08:30:14 EST 2011
Author: martinl
Date: 2011-11-11 05:30:14 -0800 (Fri, 11 Nov 2011)
New Revision: 49179
Modified:
grass/trunk/lib/vector/Vlib/build.c
grass/trunk/lib/vector/Vlib/build_nat.c
Log:
vlib: move generic fns from build_nat.c to build.c
Modified: grass/trunk/lib/vector/Vlib/build.c
===================================================================
--- grass/trunk/lib/vector/Vlib/build.c 2011-11-11 13:24:52 UTC (rev 49178)
+++ grass/trunk/lib/vector/Vlib/build.c 2011-11-11 13:30:14 UTC (rev 49179)
@@ -42,6 +42,438 @@
};
/*!
+ \brief Build area on given side of line (GV_LEFT or GV_RIGHT)
+
+ \param Map pointer to Map_info structure
+ \param iline line id
+ \param side side (GV_LEFT or GV_RIGHT)
+
+ \return > 0 number of area
+ \return < 0 number of isle
+ \return 0 not created (may also already exist)
+ */
+int Vect_build_line_area(struct Map_info *Map, int iline, int side)
+{
+ int j, area, isle, n_lines, line, direction;
+ static int first = 1;
+ off_t offset;
+ struct Plus_head *plus;
+ struct P_line *BLine;
+ static struct line_pnts *Points, *APoints;
+ struct bound_box box;
+ plus_t *lines;
+ double area_size;
+
+ plus = &(Map->plus);
+
+ G_debug(3, "Vect_build_line_area() line = %d, side = %d", iline, side);
+
+ if (first) {
+ Points = Vect_new_line_struct();
+ APoints = Vect_new_line_struct();
+ first = 0;
+ }
+
+ area = dig_line_get_area(plus, iline, side);
+ if (area != 0) {
+ G_debug(3, " area/isle = %d -> skip", area);
+ return 0;
+ }
+
+ n_lines = dig_build_area_with_line(plus, iline, side, &lines);
+ G_debug(3, " n_lines = %d", n_lines);
+ if (n_lines < 1) {
+ return 0;
+ } /* area was not built */
+
+ /* Area or island ? */
+ Vect_reset_line(APoints);
+ for (j = 0; j < n_lines; j++) {
+ line = abs(lines[j]);
+ BLine = plus->Line[line];
+ offset = BLine->offset;
+ G_debug(3, " line[%d] = %d, offset = %lu", j, line,
+ (unsigned long)offset);
+ Vect_read_line(Map, Points, NULL, line);
+ if (lines[j] > 0)
+ direction = GV_FORWARD;
+ else
+ direction = GV_BACKWARD;
+ Vect_append_points(APoints, Points, direction);
+ APoints->n_points--; /* skip last point, avoids duplicates */
+ }
+ dig_line_box(APoints, &box);
+ APoints->n_points++; /* close polygon */
+
+ dig_find_area_poly(APoints, &area_size);
+
+ /* area_size = dig_find_poly_orientation(APoints); */
+ /* area_size is not real area size, we are only interested in the sign */
+
+ G_debug(3, " area/isle size = %f", area_size);
+
+ if (area_size > 0) { /* CW: area */
+ /* add area structure to plus */
+ area = dig_add_area(plus, n_lines, lines, &box);
+ if (area == -1) { /* error */
+ Vect_close(Map);
+ G_fatal_error(_("Unable to add area (map closed, topo saved)"));
+ }
+ G_debug(3, " -> area %d", area);
+ return area;
+ }
+ else if (area_size < 0) { /* CCW: island */
+ isle = dig_add_isle(plus, n_lines, lines, &box);
+ if (isle == -1) { /* error */
+ Vect_close(Map);
+ G_fatal_error(_("Unable to add isle (map closed, topo saved)"));
+ }
+ G_debug(3, " -> isle %d", isle);
+ return -isle;
+ }
+ else {
+ /* TODO: What to do with such areas? Should be areas/isles of size 0 stored,
+ * so that may be found and cleaned by some utility
+ * Note: it would be useful for vertical closed polygons, but such would be added twice
+ * as area */
+ G_warning(_("Area of size = 0.0 ignored"));
+ }
+ return 0;
+}
+
+/*!
+ \brief Find area outside island
+
+ \param Map_info vector map
+ \param isle isle id
+
+ \return area id
+ \return 0 if not found
+ */
+int Vect_isle_find_area(struct Map_info *Map, int isle)
+{
+ int j, line, sel_area, area, poly;
+ static int first_call = 1;
+ const struct Plus_head *plus;
+ struct P_line *Line;
+ struct P_node *Node;
+ struct P_isle *Isle;
+ struct P_area *Area;
+ struct P_topo_b *topo;
+ double size, cur_size;
+ struct bound_box box, abox;
+ static struct boxlist *List;
+ static struct line_pnts *APoints;
+
+ /* Note: We should check all isle points (at least) because if topology is not clean
+ * and two areas overlap, isle which is not completely within area may be attached,
+ * but it would take long time */
+
+ G_debug(3, "Vect_isle_find_area () island = %d", isle);
+ plus = &(Map->plus);
+
+ if (plus->Isle[isle] == NULL) {
+ G_warning(_("Request to find area outside nonexistent isle"));
+ return 0;
+ }
+
+ if (first_call) {
+ List = Vect_new_boxlist(1);
+ APoints = Vect_new_line_struct();
+ first_call = 0;
+ }
+
+ Isle = plus->Isle[isle];
+ line = abs(Isle->lines[0]);
+ Line = plus->Line[line];
+ topo = (struct P_topo_b *)Line->topo;
+ Node = plus->Node[topo->N1];
+
+ /* select areas by box */
+ box.E = Node->x;
+ box.W = Node->x;
+ box.N = Node->y;
+ box.S = Node->y;
+ box.T = PORT_DOUBLE_MAX;
+ box.B = -PORT_DOUBLE_MAX;
+ Vect_select_areas_by_box(Map, &box, List);
+ G_debug(3, "%d areas overlap island boundary point", List->n_values);
+
+ sel_area = 0;
+ cur_size = -1;
+ Vect_get_isle_box(Map, isle, &box);
+ for (j = 0; j < List->n_values; j++) {
+ area = List->id[j];
+ G_debug(3, "area = %d", area);
+
+ Area = plus->Area[area];
+
+ /* Before other tests, simply exclude those areas inside isolated isles formed by one boundary */
+ if (abs(Isle->lines[0]) == abs(Area->lines[0])) {
+ G_debug(3, " area inside isolated isle");
+ continue;
+ }
+
+ /* Check box */
+ /* Note: If build is run on large files of areas imported from nontopo format (shapefile)
+ * attaching of isles takes very long time because each area is also isle and select by
+ * box all overlapping areas selects all areas with box overlapping first node.
+ * Then reading coordinates for all those areas would take a long time -> check first
+ * if isle's box is completely within area box */
+
+ abox = List->box[j];
+
+ if (box.E > abox.E || box.W < abox.W || box.N > abox.N ||
+ box.S < abox.S) {
+ G_debug(3, " isle not completely inside area box");
+ continue;
+ }
+
+ poly = Vect_point_in_area_outer_ring(Node->x, Node->y, Map, area, abox);
+ G_debug(3, " poly = %d", poly);
+
+ if (poly == 1) { /* point in area, but node is not part of area inside isle (would be poly == 2) */
+ /* In rare case island is inside more areas in that case we have to calculate area
+ * of outer ring and take the smaller */
+ if (sel_area == 0) { /* first */
+ sel_area = area;
+ }
+ else { /* is not first */
+ if (cur_size < 0) { /* second area */
+ /* This is slow, but should not be called often */
+ Vect_get_area_points(Map, sel_area, APoints);
+ /* G_begin_polygon_area_calculations();
+ cur_size =
+ G_area_of_polygon(APoints->x, APoints->y,
+ APoints->n_points); */
+ /* this is faster, but there may be latlon problems: the poles */
+ dig_find_area_poly(APoints, &cur_size);
+ G_debug(3, " first area size = %f (n points = %d)",
+ cur_size, APoints->n_points);
+
+ }
+
+ Vect_get_area_points(Map, area, APoints);
+ /* size =
+ G_area_of_polygon(APoints->x, APoints->y,
+ APoints->n_points); */
+ /* this is faster, but there may be latlon problems: the poles */
+ dig_find_area_poly(APoints, &size);
+ G_debug(3, " area size = %f (n points = %d)", size,
+ APoints->n_points);
+
+ if (size > 0 && size < cur_size) {
+ sel_area = area;
+ cur_size = size;
+ }
+ }
+ G_debug(3, "sel_area = %d cur_size = %f", sel_area, cur_size);
+ }
+ }
+ if (sel_area > 0) {
+ G_debug(3, "Island %d in area %d", isle, sel_area);
+ }
+ else {
+ G_debug(3, "Island %d is not in area", isle);
+ }
+
+ return sel_area;
+}
+
+/*!
+ \brief (Re)Attach isle to area
+
+ \param Map_info vector map
+ \param isle isle id
+
+ \return 0
+ */
+int Vect_attach_isle(struct Map_info *Map, int isle)
+{
+ int sel_area;
+ struct P_isle *Isle;
+ struct Plus_head *plus;
+
+ /* Note!: If topology is not clean and areas overlap, one island
+ may fall to more areas (partially or fully). Before isle is
+ attached to area it must be check if it is not attached yet */
+ G_debug(3, "Vect_attach_isle(): isle = %d", isle);
+
+ plus = &(Map->plus);
+
+ sel_area = Vect_isle_find_area(Map, isle);
+ G_debug(3, "\tisle = %d -> area outside = %d", isle, sel_area);
+ if (sel_area > 0) {
+ Isle = plus->Isle[isle];
+ if (Isle->area > 0) {
+ G_debug(3, "Attempt to attach isle %d to more areas "
+ "(=>topology is not clean)", isle);
+ }
+ else {
+ Isle->area = sel_area;
+ dig_area_add_isle(plus, sel_area, isle);
+ }
+ }
+ return 0;
+}
+
+/*!
+ \brief (Re)Attach isles to areas in given bounding box
+
+ \param Map_info vector map
+ \param box bounding box
+
+ \return 0
+ */
+int Vect_attach_isles(struct Map_info *Map, const struct bound_box *box)
+{
+ int i, isle;
+ static int first = TRUE;
+ static struct boxlist *List;
+ struct Plus_head *plus;
+
+ G_debug(3, "Vect_attach_isles()");
+
+ plus = &(Map->plus);
+
+ if (first) {
+ List = Vect_new_boxlist(FALSE);
+ first = FALSE;
+ }
+
+ Vect_select_isles_by_box(Map, box, List);
+ G_debug(3, " number of isles to attach = %d", List->n_values);
+
+ for (i = 0; i < List->n_values; i++) {
+ isle = List->id[i];
+ /* only attach isles that are not yet attached, see Vect_attach_isle() */
+ if (plus->Isle[isle]->area == 0)
+ Vect_attach_isle(Map, isle);
+ }
+ return 0;
+}
+
+/*!
+ \brief (Re)Attach centroids to areas in given bounding box
+
+ Warning: If map is updated on level2, it may happen that
+ previously correct island becomes incorrect. In that case,
+ centroid of area forming the island is reattached to outer area,
+ because island polygon is not excluded.
+
+ <pre>
+ +-----------+ +-----------+
+ | 1 | | 1 |
+ | +---+---+ | | +---+---+ |
+ | | 2 | 3 | | | | 2 | |
+ | | x | | | -> | | x | |
+ | | | | | | | | |
+ | +---+---+ | | +---+---+ |
+ | | | |
+ +-----------+ +-----------+
+ centroid is centroid is
+ attached to 2 reattached to 1
+ </pre>
+
+ Because of this, when the centroid is reattached to another area,
+ it is always necessary to check if original area exist, unregister
+ centroid from previous area. To simplify code, this is
+ implemented so that centroid is always firs unregistered and if
+ new area is found, it is registered again.
+
+ This problem can be avoided altogether if properly attached
+ centroids are skipped MM 2009
+
+ \param Map_info vector map
+ \param box bounding box
+
+ \return 0
+ */
+int Vect_attach_centroids(struct Map_info *Map, const struct bound_box * box)
+{
+ int i, sel_area, centr;
+ static int first = 1;
+ static struct boxlist *List;
+ static struct line_pnts *Points;
+ struct P_area *Area;
+ struct P_line *Line;
+ struct P_topo_c *topo;
+ struct Plus_head *plus;
+
+ G_debug(3, "Vect_attach_centroids()");
+
+ plus = &(Map->plus);
+
+ if (first) {
+ List = Vect_new_boxlist(0);
+ Points = Vect_new_line_struct();
+ first = 0;
+ }
+
+ Vect_select_lines_by_box(Map, box, GV_CENTROID, List);
+ G_debug(3, "\tnumber of centroids to reattach = %d", List->n_values);
+ for (i = 0; i < List->n_values; i++) {
+ int orig_area;
+
+ centr = List->id[i];
+ Line = plus->Line[centr];
+ topo = (struct P_topo_c *)Line->topo;
+
+ /* only attach unregistered and duplicate centroids because
+ * 1) all properly attached centroids are properly attached, really! Don't touch.
+ * 2) Vect_find_area() below does not always return the correct area
+ * 3) it's faster
+ */
+ if (topo->area > 0)
+ continue;
+
+ orig_area = topo->area;
+
+ Vect_read_line(Map, Points, NULL, centr);
+ if (Points->n_points < 1) {
+ /* try to get centroid from spatial index (OGR layers) */
+ int found;
+ struct boxlist list;
+ dig_init_boxlist(&list, TRUE);
+ Vect_select_lines_by_box(Map, box, GV_CENTROID, &list);
+
+ found = 0;
+ for (i = 0; i < list.n_values; i++) {
+ if (list.id[i] == centr) {
+ found = i;
+ break;
+ }
+ }
+ Vect_append_point(Points, list.box[found].E, list.box[found].N, 0.0);
+ }
+ sel_area = Vect_find_area(Map, Points->x[0], Points->y[0]);
+ G_debug(3, "\tcentroid %d is in area %d", centr, sel_area);
+ if (sel_area > 0) {
+ Area = plus->Area[sel_area];
+ if (Area->centroid == 0) { /* first centroid */
+ G_debug(3, "\tfirst centroid -> attach to area");
+ Area->centroid = centr;
+ topo->area = sel_area;
+
+ if (sel_area != orig_area && plus->do_uplist)
+ dig_line_add_updated(plus, centr);
+ }
+ else if (Area->centroid != centr) { /* duplicate centroid */
+ /* Note: it cannot happen that Area->centroid == centr, because the centroid
+ * was not registered or a duplicate */
+ G_debug(3, "\tduplicate centroid -> do not attach to area");
+ topo->area = -sel_area;
+
+ if (-sel_area != orig_area && plus->do_uplist)
+ dig_line_add_updated(plus, centr);
+ }
+ }
+ }
+
+ return 0;
+}
+
+/*!
\brief Build topology for vector map
\param Map vector map
Modified: grass/trunk/lib/vector/Vlib/build_nat.c
===================================================================
--- grass/trunk/lib/vector/Vlib/build_nat.c 2011-11-11 13:24:52 UTC (rev 49178)
+++ grass/trunk/lib/vector/Vlib/build_nat.c 2011-11-11 13:30:14 UTC (rev 49179)
@@ -3,7 +3,7 @@
\brief Vector library - Building topology for native format
- (C) 2001-2009 by the GRASS Development Team
+ (C) 2001-2009, 2011 by the GRASS Development Team
This program is free software under the
GNU General Public License (>=v2).
@@ -22,438 +22,6 @@
#include <grass/vector.h>
/*!
- \brief Build area on given side of line (GV_LEFT or GV_RIGHT)
-
- \param Map pointer to Map_info structure
- \param iline line id
- \param side side (GV_LEFT or GV_RIGHT)
-
- \return > 0 number of area
- \return < 0 number of isle
- \return 0 not created (may also already exist)
- */
-int Vect_build_line_area(struct Map_info *Map, int iline, int side)
-{
- int j, area, isle, n_lines, line, direction;
- static int first = 1;
- off_t offset;
- struct Plus_head *plus;
- struct P_line *BLine;
- static struct line_pnts *Points, *APoints;
- struct bound_box box;
- plus_t *lines;
- double area_size;
-
- plus = &(Map->plus);
-
- G_debug(3, "Vect_build_line_area() line = %d, side = %d", iline, side);
-
- if (first) {
- Points = Vect_new_line_struct();
- APoints = Vect_new_line_struct();
- first = 0;
- }
-
- area = dig_line_get_area(plus, iline, side);
- if (area != 0) {
- G_debug(3, " area/isle = %d -> skip", area);
- return 0;
- }
-
- n_lines = dig_build_area_with_line(plus, iline, side, &lines);
- G_debug(3, " n_lines = %d", n_lines);
- if (n_lines < 1) {
- return 0;
- } /* area was not built */
-
- /* Area or island ? */
- Vect_reset_line(APoints);
- for (j = 0; j < n_lines; j++) {
- line = abs(lines[j]);
- BLine = plus->Line[line];
- offset = BLine->offset;
- G_debug(3, " line[%d] = %d, offset = %lu", j, line,
- (unsigned long)offset);
- Vect_read_line(Map, Points, NULL, line);
- if (lines[j] > 0)
- direction = GV_FORWARD;
- else
- direction = GV_BACKWARD;
- Vect_append_points(APoints, Points, direction);
- APoints->n_points--; /* skip last point, avoids duplicates */
- }
- dig_line_box(APoints, &box);
- APoints->n_points++; /* close polygon */
-
- dig_find_area_poly(APoints, &area_size);
-
- /* area_size = dig_find_poly_orientation(APoints); */
- /* area_size is not real area size, we are only interested in the sign */
-
- G_debug(3, " area/isle size = %f", area_size);
-
- if (area_size > 0) { /* CW: area */
- /* add area structure to plus */
- area = dig_add_area(plus, n_lines, lines, &box);
- if (area == -1) { /* error */
- Vect_close(Map);
- G_fatal_error(_("Unable to add area (map closed, topo saved)"));
- }
- G_debug(3, " -> area %d", area);
- return area;
- }
- else if (area_size < 0) { /* CCW: island */
- isle = dig_add_isle(plus, n_lines, lines, &box);
- if (isle == -1) { /* error */
- Vect_close(Map);
- G_fatal_error(_("Unable to add isle (map closed, topo saved)"));
- }
- G_debug(3, " -> isle %d", isle);
- return -isle;
- }
- else {
- /* TODO: What to do with such areas? Should be areas/isles of size 0 stored,
- * so that may be found and cleaned by some utility
- * Note: it would be useful for vertical closed polygons, but such would be added twice
- * as area */
- G_warning(_("Area of size = 0.0 ignored"));
- }
- return 0;
-}
-
-/*!
- \brief Find area outside island
-
- \param Map_info vector map
- \param isle isle id
-
- \return area id
- \return 0 if not found
- */
-int Vect_isle_find_area(struct Map_info *Map, int isle)
-{
- int j, line, sel_area, area, poly;
- static int first_call = 1;
- const struct Plus_head *plus;
- struct P_line *Line;
- struct P_node *Node;
- struct P_isle *Isle;
- struct P_area *Area;
- struct P_topo_b *topo;
- double size, cur_size;
- struct bound_box box, abox;
- static struct boxlist *List;
- static struct line_pnts *APoints;
-
- /* Note: We should check all isle points (at least) because if topology is not clean
- * and two areas overlap, isle which is not completely within area may be attached,
- * but it would take long time */
-
- G_debug(3, "Vect_isle_find_area () island = %d", isle);
- plus = &(Map->plus);
-
- if (plus->Isle[isle] == NULL) {
- G_warning(_("Request to find area outside nonexistent isle"));
- return 0;
- }
-
- if (first_call) {
- List = Vect_new_boxlist(1);
- APoints = Vect_new_line_struct();
- first_call = 0;
- }
-
- Isle = plus->Isle[isle];
- line = abs(Isle->lines[0]);
- Line = plus->Line[line];
- topo = (struct P_topo_b *)Line->topo;
- Node = plus->Node[topo->N1];
-
- /* select areas by box */
- box.E = Node->x;
- box.W = Node->x;
- box.N = Node->y;
- box.S = Node->y;
- box.T = PORT_DOUBLE_MAX;
- box.B = -PORT_DOUBLE_MAX;
- Vect_select_areas_by_box(Map, &box, List);
- G_debug(3, "%d areas overlap island boundary point", List->n_values);
-
- sel_area = 0;
- cur_size = -1;
- Vect_get_isle_box(Map, isle, &box);
- for (j = 0; j < List->n_values; j++) {
- area = List->id[j];
- G_debug(3, "area = %d", area);
-
- Area = plus->Area[area];
-
- /* Before other tests, simply exclude those areas inside isolated isles formed by one boundary */
- if (abs(Isle->lines[0]) == abs(Area->lines[0])) {
- G_debug(3, " area inside isolated isle");
- continue;
- }
-
- /* Check box */
- /* Note: If build is run on large files of areas imported from nontopo format (shapefile)
- * attaching of isles takes very long time because each area is also isle and select by
- * box all overlapping areas selects all areas with box overlapping first node.
- * Then reading coordinates for all those areas would take a long time -> check first
- * if isle's box is completely within area box */
-
- abox = List->box[j];
-
- if (box.E > abox.E || box.W < abox.W || box.N > abox.N ||
- box.S < abox.S) {
- G_debug(3, " isle not completely inside area box");
- continue;
- }
-
- poly = Vect_point_in_area_outer_ring(Node->x, Node->y, Map, area, abox);
- G_debug(3, " poly = %d", poly);
-
- if (poly == 1) { /* point in area, but node is not part of area inside isle (would be poly == 2) */
- /* In rare case island is inside more areas in that case we have to calculate area
- * of outer ring and take the smaller */
- if (sel_area == 0) { /* first */
- sel_area = area;
- }
- else { /* is not first */
- if (cur_size < 0) { /* second area */
- /* This is slow, but should not be called often */
- Vect_get_area_points(Map, sel_area, APoints);
- /* G_begin_polygon_area_calculations();
- cur_size =
- G_area_of_polygon(APoints->x, APoints->y,
- APoints->n_points); */
- /* this is faster, but there may be latlon problems: the poles */
- dig_find_area_poly(APoints, &cur_size);
- G_debug(3, " first area size = %f (n points = %d)",
- cur_size, APoints->n_points);
-
- }
-
- Vect_get_area_points(Map, area, APoints);
- /* size =
- G_area_of_polygon(APoints->x, APoints->y,
- APoints->n_points); */
- /* this is faster, but there may be latlon problems: the poles */
- dig_find_area_poly(APoints, &size);
- G_debug(3, " area size = %f (n points = %d)", size,
- APoints->n_points);
-
- if (size > 0 && size < cur_size) {
- sel_area = area;
- cur_size = size;
- }
- }
- G_debug(3, "sel_area = %d cur_size = %f", sel_area, cur_size);
- }
- }
- if (sel_area > 0) {
- G_debug(3, "Island %d in area %d", isle, sel_area);
- }
- else {
- G_debug(3, "Island %d is not in area", isle);
- }
-
- return sel_area;
-}
-
-/*!
- \brief (Re)Attach isle to area
-
- \param Map_info vector map
- \param isle isle id
-
- \return 0
- */
-int Vect_attach_isle(struct Map_info *Map, int isle)
-{
- int sel_area;
- struct P_isle *Isle;
- struct Plus_head *plus;
-
- /* Note!: If topology is not clean and areas overlap, one island
- may fall to more areas (partially or fully). Before isle is
- attached to area it must be check if it is not attached yet */
- G_debug(3, "Vect_attach_isle(): isle = %d", isle);
-
- plus = &(Map->plus);
-
- sel_area = Vect_isle_find_area(Map, isle);
- G_debug(3, "\tisle = %d -> area outside = %d", isle, sel_area);
- if (sel_area > 0) {
- Isle = plus->Isle[isle];
- if (Isle->area > 0) {
- G_debug(3, "Attempt to attach isle %d to more areas "
- "(=>topology is not clean)", isle);
- }
- else {
- Isle->area = sel_area;
- dig_area_add_isle(plus, sel_area, isle);
- }
- }
- return 0;
-}
-
-/*!
- \brief (Re)Attach isles to areas in given bounding box
-
- \param Map_info vector map
- \param box bounding box
-
- \return 0
- */
-int Vect_attach_isles(struct Map_info *Map, const struct bound_box *box)
-{
- int i, isle;
- static int first = TRUE;
- static struct boxlist *List;
- struct Plus_head *plus;
-
- G_debug(3, "Vect_attach_isles()");
-
- plus = &(Map->plus);
-
- if (first) {
- List = Vect_new_boxlist(FALSE);
- first = FALSE;
- }
-
- Vect_select_isles_by_box(Map, box, List);
- G_debug(3, " number of isles to attach = %d", List->n_values);
-
- for (i = 0; i < List->n_values; i++) {
- isle = List->id[i];
- /* only attach isles that are not yet attached, see Vect_attach_isle() */
- if (plus->Isle[isle]->area == 0)
- Vect_attach_isle(Map, isle);
- }
- return 0;
-}
-
-/*!
- \brief (Re)Attach centroids to areas in given bounding box
-
- Warning: If map is updated on level2, it may happen that
- previously correct island becomes incorrect. In that case,
- centroid of area forming the island is reattached to outer area,
- because island polygon is not excluded.
-
- <pre>
- +-----------+ +-----------+
- | 1 | | 1 |
- | +---+---+ | | +---+---+ |
- | | 2 | 3 | | | | 2 | |
- | | x | | | -> | | x | |
- | | | | | | | | |
- | +---+---+ | | +---+---+ |
- | | | |
- +-----------+ +-----------+
- centroid is centroid is
- attached to 2 reattached to 1
- </pre>
-
- Because of this, when the centroid is reattached to another area,
- it is always necessary to check if original area exist, unregister
- centroid from previous area. To simplify code, this is
- implemented so that centroid is always firs unregistered and if
- new area is found, it is registered again.
-
- This problem can be avoided altogether if properly attached
- centroids are skipped MM 2009
-
- \param Map_info vector map
- \param box bounding box
-
- \return 0
- */
-int Vect_attach_centroids(struct Map_info *Map, const struct bound_box * box)
-{
- int i, sel_area, centr;
- static int first = 1;
- static struct boxlist *List;
- static struct line_pnts *Points;
- struct P_area *Area;
- struct P_line *Line;
- struct P_topo_c *topo;
- struct Plus_head *plus;
-
- G_debug(3, "Vect_attach_centroids()");
-
- plus = &(Map->plus);
-
- if (first) {
- List = Vect_new_boxlist(0);
- Points = Vect_new_line_struct();
- first = 0;
- }
-
- Vect_select_lines_by_box(Map, box, GV_CENTROID, List);
- G_debug(3, "\tnumber of centroids to reattach = %d", List->n_values);
- for (i = 0; i < List->n_values; i++) {
- int orig_area;
-
- centr = List->id[i];
- Line = plus->Line[centr];
- topo = (struct P_topo_c *)Line->topo;
-
- /* only attach unregistered and duplicate centroids because
- * 1) all properly attached centroids are properly attached, really! Don't touch.
- * 2) Vect_find_area() below does not always return the correct area
- * 3) it's faster
- */
- if (topo->area > 0)
- continue;
-
- orig_area = topo->area;
-
- Vect_read_line(Map, Points, NULL, centr);
- if (Points->n_points < 1) {
- /* try to get centroid from spatial index (OGR layers) */
- int found;
- struct boxlist list;
- dig_init_boxlist(&list, TRUE);
- Vect_select_lines_by_box(Map, box, GV_CENTROID, &list);
-
- found = 0;
- for (i = 0; i < list.n_values; i++) {
- if (list.id[i] == centr) {
- found = i;
- break;
- }
- }
- Vect_append_point(Points, list.box[found].E, list.box[found].N, 0.0);
- }
- sel_area = Vect_find_area(Map, Points->x[0], Points->y[0]);
- G_debug(3, "\tcentroid %d is in area %d", centr, sel_area);
- if (sel_area > 0) {
- Area = plus->Area[sel_area];
- if (Area->centroid == 0) { /* first centroid */
- G_debug(3, "\tfirst centroid -> attach to area");
- Area->centroid = centr;
- topo->area = sel_area;
-
- if (sel_area != orig_area && plus->do_uplist)
- dig_line_add_updated(plus, centr);
- }
- else if (Area->centroid != centr) { /* duplicate centroid */
- /* Note: it cannot happen that Area->centroid == centr, because the centroid
- * was not registered or a duplicate */
- G_debug(3, "\tduplicate centroid -> do not attach to area");
- topo->area = -sel_area;
-
- if (-sel_area != orig_area && plus->do_uplist)
- dig_line_add_updated(plus, centr);
- }
- }
- }
-
- return 0;
-}
-
-/*!
\brief Build topology
\param Map_info vector map
More information about the grass-commit
mailing list