vlib: move generic fns from build_nat.c to build.c

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@49179 15284696-431f-4ddb-bdfa-cd5b030d7da7

lib/vector/Vlib/build.c

@@ -42,6 +42,438 @@ static int (*Build_array[]) () = {
 };
 
 /*!
+   \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

lib/vector/Vlib/build_nat.c

@@ -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