only attach isles that are not attached

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

lib/vector/Vlib/build_nat.c

@@ -324,7 +324,9 @@ int Vect_attach_isles(struct Map_info *Map, const struct bound_box * box)
 
     for (i = 0; i < List->n_values; i++) {
 	isle = List->value[i];
-	Vect_attach_isle(Map, isle);
+	/* 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;
 }

lib/vector/Vlib/remove_areas.c

@@ -104,7 +104,7 @@ Vect_remove_small_areas(struct Map_info *Map, double thresh,
 	}
 	G_debug(3, "num neighbours = %d", AList->n_values);
 
-	/* Go through the list of neghours and find that with the longest boundary */
+	/* Go through the list of neighours and find that with the longest boundary */
 	dissolve_neighbour = 0;
 	length = -1.0;
 	for (i = 0; i < AList->n_values; i++) {
@@ -178,5 +178,483 @@ Vect_remove_small_areas(struct Map_info *Map, double thresh,
     G_verbose_message(_("%d areas of total size %g removed"), nremoved,
 		size_removed);
 
+    Vect_build_partial(Map, GV_BUILD_BASE);
+    Vect_merge_lines(Map, GV_BOUNDARY, NULL, NULL);
+
+    return (nremoved);
+}
+
+/* faster version */
+int
+Vect_remove_small_areas_faster(struct Map_info *Map, double thresh,
+			struct Map_info *Err, double *removed_area)
+{
+    int area, nareas;
+    int nremoved = 0;
+    struct ilist *List;
+    struct ilist *AList;
+    struct line_pnts *Points;
+    struct line_cats *Cats;
+    double size_removed = 0.0;
+    int left, right, dissolve_neighbour;
+
+    List = Vect_new_list();
+    AList = Vect_new_list();
+    Points = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+
+    nareas = Vect_get_num_areas(Map);
+    for (area = 1; area <= Vect_get_num_areas(Map); area++) {
+	int i, centroid, remove_boundary;
+	double length, size;
+
+	if (area <= nareas)
+	    G_percent(area, nareas, 1);
+	G_debug(3, "area = %d", area);
+	if (!Vect_area_alive(Map, area))
+	    continue;
+
+	size = Vect_get_area_area(Map, area);
+	if (size > thresh)
+	    continue;
+	size_removed += size;
+
+	/* The area is smaller than the limit -> remove */
+
+	/* Remove centroid */
+	centroid = Vect_get_area_centroid(Map, area);
+	if (centroid > 0) {
+	    if (Err) {
+		Vect_read_line(Map, Points, Cats, centroid);
+		Vect_write_line(Err, GV_CENTROID, Points, Cats);
+	    }
+	    Vect_delete_line(Map, centroid);
+	}
+
+	/* Find the adjacent area with which the longest boundary is shared */
+
+	Vect_get_area_boundaries(Map, area, List);
+
+	/* Go through the list of boundaries and find the longest boundary */
+	remove_boundary = 0;
+	length = -1.0;
+	for (i = 0; i < List->n_values; i++) {
+	    int line;
+	    double l = 0.0;
+
+	    line = List->value[i];
+	    G_debug(4, "   line = %d", line);
+
+	    Vect_read_line(Map, Points, NULL, abs(line));
+	    l = Vect_line_length(Points);
+
+	    if (l > length) {
+		length = l;
+		remove_boundary = line;
+	    }
+	}
+
+	G_debug(3, "remove_boundary = %d", remove_boundary);
+
+	Vect_get_line_areas(Map, abs(remove_boundary), &left, &right);
+	dissolve_neighbour = 0;
+	if (remove_boundary > 0)
+	    dissolve_neighbour = left;
+	else
+	    dissolve_neighbour = right;
+
+	G_debug(3, "dissolve_neighbour = %d", dissolve_neighbour);
+	G_debug(3, "dissolve_neighbour = %d", dissolve_neighbour);
+	if (dissolve_neighbour == 0)
+	    G_fatal_error("could not find neighbour to dissolve");
+
+	/* Make list of boundaries to be removed */
+	Vect_reset_list(AList);
+	for (i = 0; i < List->n_values; i++) {
+	    int line, neighbour;
+
+	    line = List->value[i];
+	    Vect_get_line_areas(Map, abs(line), &left, &right);
+	    if (line > 0)
+		neighbour = left;
+	    else
+		neighbour = right;
+
+	    G_debug(3, "   neighbour = %d", neighbour);
+
+	    if (neighbour == dissolve_neighbour) {
+		Vect_list_append(AList, abs(line));
+	    }
+	}
+
+	/* Remove boundaries */
+	for (i = 0; i < AList->n_values; i++) {
+	    int line;
+
+	    line = AList->value[i];
+
+	    if (Err) {
+		Vect_read_line(Map, Points, Cats, line);
+		Vect_write_line(Err, GV_BOUNDARY, Points, Cats);
+	    }
+	    Vect_delete_line(Map, line);
+	}
+
+	nremoved++;
+    }
+
+    if (removed_area)
+	*removed_area = size_removed;
+
+    G_verbose_message(_("%d areas of total size %g removed"), nremoved,
+		size_removed);
+
+    Vect_build_partial(Map, GV_BUILD_BASE);
+    Vect_merge_lines(Map, GV_BOUNDARY, NULL, NULL);
+
+    return (nremoved);
+}
+
+/* much faster version */
+int
+Vect_remove_small_areas_fastest(struct Map_info *Map, double thresh,
+			struct Map_info *Err, double *removed_area)
+{
+    int area, nareas;
+    int nremoved = 0;
+    struct ilist *List;
+    struct ilist *AList;
+    struct ilist *BList;
+    struct ilist *NList;
+    struct ilist *IList;
+    struct line_pnts *Points;
+    struct line_cats *Cats;
+    double size_removed = 0.0;
+    int dissolve_neighbour;
+    int line, left, right, neighbour;
+    int nisles, nnisles;
+
+    List = Vect_new_list();
+    AList = Vect_new_list();
+    BList = Vect_new_list();
+    NList = Vect_new_list();
+    IList = Vect_new_list();
+    Points = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+
+    nareas = Vect_get_num_areas(Map);
+    for (area = 1; area <= Vect_get_num_areas(Map); area++) {
+	int i, j, centroid;
+	double length, l, size;
+	int outer_isle = 1, outer_area = -1;
+
+	if (area <= nareas)
+	    G_percent(area, nareas, 1);
+	G_debug(3, "area = %d", area);
+	if (!Vect_area_alive(Map, area))
+	    continue;
+
+	size = Vect_get_area_area(Map, area);
+	if (size > thresh)
+	    continue;
+	size_removed += size;
+
+	/* The area is smaller than the limit -> remove */
+
+	/* Remove centroid */
+	centroid = Vect_get_area_centroid(Map, area);
+	if (centroid > 0) {
+	    if (Err) {
+		Vect_read_line(Map, Points, Cats, centroid);
+		Vect_write_line(Err, GV_CENTROID, Points, Cats);
+	    }
+	    Vect_delete_line(Map, centroid);
+	}
+
+	/* Find the adjacent area with which the longest boundary is shared */
+
+	Vect_get_area_boundaries(Map, area, List);
+
+	/* Create a list of neighbour areas */
+	Vect_reset_list(AList);
+	for (i = 0; i < List->n_values; i++) {
+
+	    line = List->value[i];
+
+	    if (!Vect_line_alive(Map, abs(line)))	/* Should not happen */
+		G_fatal_error(_("Area is composed of dead boundary"));
+
+	    Vect_get_line_areas(Map, abs(line), &left, &right);
+	    if (line > 0)
+		neighbour = left;
+	    else
+		neighbour = right;
+
+	    G_debug(4, "  line = %d left = %d right = %d neighbour = %d",
+		    line, left, right, neighbour);
+
+	    Vect_list_append(AList, neighbour);	/* this checks for duplicity */
+	}
+	G_debug(3, "num neighbours = %d", AList->n_values);
+
+	/* Go through the list of neighbours and find the one with the longest boundary */
+	dissolve_neighbour = 0;
+	length = -1.0;
+	for (i = 0; i < AList->n_values; i++) {
+	    int neighbour1;
+
+	    l = 0.0;
+	    neighbour1 = AList->value[i];
+	    G_debug(4, "   neighbour1 = %d", neighbour1);
+
+	    for (j = 0; j < List->n_values; j++) {
+		int neighbour2;
+
+		line = List->value[j];
+		Vect_get_line_areas(Map, abs(line), &left, &right);
+		if (line > 0)
+		    neighbour2 = left;
+		else
+		    neighbour2 = right;
+
+		if (neighbour2 == neighbour1) {
+		    Vect_read_line(Map, Points, NULL, abs(line));
+		    l += Vect_line_length(Points);
+		}
+	    }
+	    if (l > length) {
+		length = l;
+		dissolve_neighbour = neighbour1;
+	    }
+	}
+
+	G_debug(3, "dissolve_neighbour = %d", dissolve_neighbour);
+
+	if (dissolve_neighbour == 0) {
+	    G_fatal_error("could not find neighbour to dissolve");
+	}
+
+	/* Make list of boundaries to be removed */
+	Vect_reset_list(AList);
+	Vect_reset_list(BList);
+	for (i = 0; i < List->n_values; i++) {
+
+	    line = List->value[i];
+	    Vect_get_line_areas(Map, abs(line), &left, &right);
+	    if (line > 0)
+		neighbour = left;
+	    else
+		neighbour = right;
+
+	    G_debug(3, "   neighbour = %d", neighbour);
+
+	    if (neighbour == dissolve_neighbour) {
+		Vect_list_append(AList, abs(line));
+	    }
+	    else
+		Vect_list_append(BList, line);
+		
+	    if (neighbour < 0) {
+		if (outer_isle > 0)
+		    outer_isle = neighbour;
+		else if (outer_isle != neighbour)
+		    G_fatal_error("Two different isles outside area");
+	    }
+	}
+	G_debug(3, "remove %d of %d boundaries", AList->n_values, List->n_values);
+
+	/* Get isles inside area */
+	Vect_reset_list(IList);
+	if ((nisles = Vect_get_area_num_isles(Map, area)) > 0) {
+	    for (i = 0; i < nisles; i++) {
+		Vect_list_append(IList, Vect_get_area_isle(Map, area, i));
+	    }
+	}
+
+	/* Remove boundaries */
+	for (i = 0; i < AList->n_values; i++) {
+	    int ret;
+
+	    line = AList->value[i];
+
+	    if (Err) {
+		Vect_read_line(Map, Points, Cats, line);
+		Vect_write_line(Err, GV_BOUNDARY, Points, Cats);
+	    }
+	    /* Vect_delete_line(Map, line); */
+
+	    /* delete the line from coor */
+	    ret = V1_delete_line_nat(Map, Map->plus.Line[line]->offset);
+
+	    if (ret == -1) {
+		G_fatal_error("Could not delete line from coor");
+	    }
+	}
+
+	/* update topo */
+	if (outer_isle < 0) {
+	    /* outer area */
+	    outer_area = Vect_get_isle_area(Map, -outer_isle);
+	}
+
+	if (dissolve_neighbour > 0) {
+	    int nareas_old;
+
+	    G_debug(3, "dissolve with neighbour area");
+
+	    /* get neighbour centroid */
+	    centroid = Vect_get_area_centroid(Map, dissolve_neighbour);
+	    /* get neighbour isles */
+	    if ((nnisles = Vect_get_area_num_isles(Map, dissolve_neighbour)) > 0) {
+		for (i = 0; i < nnisles; i++) {
+		    Vect_list_append(IList, Vect_get_area_isle(Map, dissolve_neighbour, i));
+		}
+	    }
+	    /* get neighbour boundaries */
+	    Vect_get_area_boundaries(Map, dissolve_neighbour, NList);
+
+	    /* delete area from topo */
+	    dig_del_area(&(Map->plus), area);
+	    /* delete neighbour area from topo */
+	    dig_del_area(&(Map->plus), dissolve_neighbour);
+	    /* delete boundaries from topo */
+	    for (i = 0; i < AList->n_values; i++) {
+		line = AList->value[i];
+		dig_del_line(&(Map->plus), line);
+	    }
+	    /* rebuild neighbour area from leftover boundaries */
+	    nareas_old = Vect_get_num_areas(Map);
+
+	    for (i = 0; i < BList->n_values; i++) {
+		line = BList->value[i];
+		if (Map->plus.Line[abs(line)]->left == 0 || Map->plus.Line[abs(line)]->right == 0) {
+		    int new_isle;
+
+		    new_isle = Vect_build_line_area(Map, abs(line), (line > 0 ? GV_RIGHT : GV_LEFT));
+		    if (new_isle > 0) {
+			outer_area = new_isle;
+			/* reattach centroid */
+			Map->plus.Area[outer_area]->centroid = centroid;
+			if (centroid > 0)
+			    Map->plus.Line[centroid]->left = outer_area;
+		    }
+		    else if (new_isle < 0) {
+			Vect_list_append(IList, -new_isle);
+		    }
+		}
+		/* check */
+		if (Map->plus.Line[abs(line)]->left == 0 || Map->plus.Line[abs(line)]->right == 0)
+		    G_fatal_error("dissolve with neighbour area: corrupt topology");
+	    }
+	    for (i = 0; i < NList->n_values; i++) {
+		line = NList->value[i];
+		if (!Vect_line_alive(Map, abs(line)))
+		    continue;
+
+		if (Map->plus.Line[abs(line)]->left == 0 || Map->plus.Line[abs(line)]->right == 0) {
+		    int new_isle;
+
+		    new_isle = Vect_build_line_area(Map, abs(line), (line > 0 ? GV_RIGHT : GV_LEFT));
+		    if (new_isle > 0) {
+			if (outer_area < 0) {
+			    outer_area = new_isle;
+			    /* reattach centroid */
+			    Map->plus.Area[outer_area]->centroid = centroid;
+			    if (centroid > 0)
+				Map->plus.Line[centroid]->left = outer_area;
+			}
+			else
+			    G_fatal_error("WTF");
+		    }
+		    else if (new_isle < 0) {
+			Vect_list_append(IList, -new_isle);
+		    }
+		}
+		if (Map->plus.Line[abs(line)]->left == 0 || Map->plus.Line[abs(line)]->right == 0)
+		    G_fatal_error("dissolve with neighbour area n bounds: corrupt topology");
+	    }
+	}
+	/* dissolve with outer isle */
+	else if (dissolve_neighbour < 0) {
+	    int nisles_old;
+
+	    G_debug(3, "dissolve with outer isle");
+
+	    if (dissolve_neighbour != outer_isle)
+		G_fatal_error("wrong outer isle to dissolve");
+
+	    /* get isle boundaries */
+	    Vect_get_isle_boundaries(Map, -dissolve_neighbour, NList);
+	    /* delete area from topo */
+	    dig_del_area(&(Map->plus), area);
+	    /* delete isle from topo */
+	    dig_del_isle(&(Map->plus), -dissolve_neighbour);
+	    /* delete boundaries from topo */
+	    for (i = 0; i < AList->n_values; i++) {
+		line = AList->value[i];
+		dig_del_line(&(Map->plus), line);
+	    }
+	    /* rebuild isles from leftover boundaries */
+	    for (i = 0; i < BList->n_values; i++) {
+		line = BList->value[i];
+
+		if (Map->plus.Line[abs(line)]->left == 0 || Map->plus.Line[abs(line)]->right == 0) {
+		    int new_isle;
+
+		    nisles_old = Vect_get_num_islands(Map);
+		    new_isle = Vect_build_line_area(Map, abs(line), (line > 0 ? GV_RIGHT : GV_LEFT));
+		    if (new_isle > 0 || (new_isle < 0 && nisles_old == Vect_get_num_islands(Map)))
+			G_fatal_error("failed to build new isle");
+
+		    if (new_isle < 0) {
+			Vect_list_append(IList, -new_isle);
+		    }
+		}
+		/* check */
+		if (Map->plus.Line[abs(line)]->left == 0 || Map->plus.Line[abs(line)]->right == 0)
+		    G_fatal_error("dissolve with outer isle: corrupt topology");
+	    }
+	    for (i = 0; i < NList->n_values; i++) {
+		line = NList->value[i];
+		if (!Vect_line_alive(Map, abs(line)))
+		    continue;
+
+		if (Map->plus.Line[abs(line)]->left == 0 || Map->plus.Line[abs(line)]->right == 0) {
+		    int new_isle;
+
+		    new_isle = Vect_build_line_area(Map, abs(line), (line > 0 ? GV_RIGHT : GV_LEFT));
+		    if (new_isle >= 0)
+			G_fatal_error("WTF");
+		    else if (new_isle < 0) {
+			Vect_list_append(IList, -new_isle);
+		    }
+		}
+		/* check */
+		if (Map->plus.Line[abs(line)]->left == 0 || Map->plus.Line[abs(line)]->right == 0)
+		    G_fatal_error("dissolve with outer isle n bounds: corrupt topology");
+	    }
+	}
+
+	/* attach all isles to outer or new area */
+	if (outer_area >= 0) {
+	    for (i = 0; i < IList->n_values; i++) {
+		Map->plus.Isle[IList->value[i]]->area = outer_area;
+		if (outer_area > 0)
+		    dig_area_add_isle(&(Map->plus), outer_area, IList->value[i]);
+	    }
+	}
+
+	nremoved++;
+    }
+
+    if (removed_area)
+	*removed_area = size_removed;
+
+    G_verbose_message(_("%d areas of total size %g removed"), nremoved,
+		size_removed);
+
+    Vect_build_partial(Map, GV_BUILD_BASE);
+    Vect_merge_lines(Map, GV_BOUNDARY, NULL, NULL);
+
     return (nremoved);
 }

lib/vector/Vlib/write_nat.c

@@ -251,7 +251,7 @@ static void V2__add_line_to_topo_nat(struct Map_info *Map, int line,
 	    /* Reattach all centroids/isles in deleted areas + new area.
 	     *  Because isles are selected by box it covers also possible new isle created above */
 	    if (!first) {	/* i.e. old area/isle was deleted or new one created */
-		/* Reattache isles */
+		/* Reattach isles */
 		if (plus->built >= GV_BUILD_ATTACH_ISLES)
 		    Vect_attach_isles(Map, &abox);
 
@@ -782,16 +782,12 @@ int V2_delete_line_nat(struct Map_info *Map, off_t line)
 
     /* Rebuild areas/isles and attach centroids and isles */
     if (plus->built >= GV_BUILD_AREAS && type == GV_BOUNDARY) {
-	int *new_areas, nnew_areas;
+	int new_areas[4], nnew_areas;
 
 	nnew_areas = 0;
-	new_areas = (int *)G_malloc(2 * n_adjacent * sizeof(int));
 	/* Rebuild areas/isles */
 	for (i = 0; i < n_adjacent; i++) {
-	    if (adjacent[i] > 0)
-		side = GV_RIGHT;
-	    else
-		side = GV_LEFT;
+	    side = adjacent[i] > 0 ? GV_RIGHT : GV_LEFT;
 
 	    G_debug(3, "Build area for line = %d, side = %d", adjacent[i],
 		    side);