add new fns for topology check

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

lib/vector/Vlib/break_lines.c

@@ -17,6 +17,10 @@
 #include <grass/vector.h>
 #include <grass/glocale.h>
 
+static int break_lines(struct Map_info *, struct ilist *, struct ilist *,
+		       int, struct Map_info *, int);
+
+
 /*!
    \brief Break lines in vector map at each intersection.
 
@@ -32,7 +36,7 @@
 void
 Vect_break_lines(struct Map_info *Map, int type, struct Map_info *Err)
 {
-    Vect_break_lines_list(Map, NULL, NULL, type, Err);
+    break_lines(Map, NULL, NULL, type, Err, 0);
 
     return;
 }
@@ -61,7 +65,60 @@ Vect_break_lines(struct Map_info *Map, int type, struct Map_info *Err)
 
 int
 Vect_break_lines_list(struct Map_info *Map, struct ilist *List_break,
-		      struct ilist *List_ref, int type, struct Map_info *Err)
+		      struct ilist *List_ref, int type,
+		      struct Map_info *Err)
+{
+    return break_lines(Map, List_break, List_ref, type, Err, 0);
+}
+
+/*!
+   \brief Check for and count intersecting lines, do not break.
+
+   For details see Vect_break_lines_list().
+
+   \param Map input vector map 
+   \param type feature type
+   \param[out] Err vector map where points at intersections will be written or NULL
+
+   \return
+ */
+
+int
+Vect_check_line_breaks(struct Map_info *Map, int type, struct Map_info *Err)
+{
+    return break_lines(Map, NULL, NULL, type, Err, 1);
+}
+
+/*!
+   \brief Check for and count intersecting lines, do not break.
+
+   If <i>List_break</i> is given, only lines in the list are checked for
+   intersections.
+
+   If reference lines are given (<i>List_ref</i>) break only lines
+   which intersect reference lines.
+
+   \param Map input vector map 
+   \param List_break list of lines (NULL for all lines in vector map)
+   \param List_ref list of reference lines or NULL
+   \param type feature type
+   \param[out] Err vector map where points at intersections will be written or NULL
+
+   \return number of intersections
+ */
+
+int
+Vect_check_line_breaks_list(struct Map_info *Map, struct ilist *List_break,
+		      struct ilist *List_ref, int type,
+		      struct Map_info *Err)
+{
+    return break_lines(Map, List_break, List_ref, type, Err, 1);
+}
+
+int
+break_lines(struct Map_info *Map, struct ilist *List_break,
+		      struct ilist *List_ref, int type,
+		      struct Map_info *Err, int check)
 {
     struct line_pnts *APoints, *BPoints, *Points;
     struct line_pnts **AXLines, **BXLines;
@@ -177,7 +234,7 @@ Vect_break_lines_list(struct Map_info *Map, struct ilist *List_break,
 
 	    btype = Vect_read_line(Map, BPoints, BCats, bline);
 
-	    /* Check if thouch by end node only */
+	    /* Check if touch by end node only */
 	    if (!is3d) {
 		Vect_get_line_nodes(Map, aline, &anode1, &anode2);
 		Vect_get_line_nodes(Map, bline, &bnode1, &bnode2);
@@ -267,13 +324,16 @@ Vect_break_lines_list(struct Map_info *Map, struct ilist *List_break,
 	    }
 	    nx = 0;		/* number of intersections to be written to Err */
 	    if (naxlines > 0) {	/* intersection -> write out */
-		Vect_delete_line(Map, aline);
+		if (!check)
+		    Vect_delete_line(Map, aline);
 		for (k = 0; k < naxlines; k++) {
 		    /* Write new line segments */
 		    /* line may collapse, don't write zero length lines */
 		    Vect_line_prune(AXLines[k]);
 		    if ((atype & GV_POINTS) || AXLines[k]->n_points > 1) {
-			ret = Vect_write_line(Map, atype, AXLines[k], ACats);
+			if (!check)
+			    ret = Vect_write_line(Map, atype, AXLines[k],
+			                          ACats);
 			if (List_ref) {
 			    Vect_list_append(List_ref, ret);
 			}
@@ -302,15 +362,17 @@ Vect_break_lines_list(struct Map_info *Map, struct ilist *List_break,
 
 	    if (nbxlines > 0) {
 		if (aline != bline) {	/* Self intersection, do not write twice, TODO: is it OK? */
-		    Vect_delete_line(Map, bline);
+		    if (!check)
+			Vect_delete_line(Map, bline);
 		    for (k = 0; k < nbxlines; k++) {
 			/* Write new line segments */
 			/* line may collapse, don't write zero length lines */
 			Vect_line_prune(BXLines[k]);
 			if ((btype & GV_POINTS) || BXLines[k]->n_points > 1) {
-			    ret =
-				Vect_write_line(Map, btype, BXLines[k],
-						BCats);
+			    if (!check)
+				ret =
+				    Vect_write_line(Map, btype, BXLines[k],
+						    BCats);
 			    G_debug(5, "Line %d written", ret);
 			    if (List_break) {
 				Vect_list_append(List_break, ret);

lib/vector/Vlib/build.c

@@ -487,6 +487,158 @@ int Vect_build(struct Map_info *Map)
 }
 
 /*!
+   \brief Extensive tests for correct topology
+
+   \param Map vector map
+   \param[out] Err vector map where errors will be written
+
+   \return 1 on success
+   \return 0 on error
+ */
+int Vect_topo_check(struct Map_info *Map, struct Map_info *Err)
+{
+    int line, nlines;
+    int nerrors, n_zero_lines, n_zero_boundaries;
+    struct line_pnts *Points;
+    struct line_cats *Cats;
+
+    if (Vect_get_built(Map) != GV_BUILD_ALL) {
+	Vect_build_partial(Map, GV_BUILD_NONE);
+	Vect_build(Map);
+    }
+
+    G_message(_("Checking for topological errors..."));
+
+    Points = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+
+    /* boundaries of zero length */
+    n_zero_lines = n_zero_boundaries = 0;
+    nlines = Vect_get_num_lines(Map);
+    for (line = 1; line <= nlines; line++) {
+	int type;
+	
+	if (!Vect_line_alive(Map, line))
+	    continue;
+	    
+	type = Vect_get_line_type(Map, line);
+
+	if (type & GV_LINES) {
+	    double len;
+	    
+	    Vect_read_line(Map, Points, Cats, line);
+	    len = Vect_line_length(Points);
+	    
+	    if (len == 0) {
+		if (type & GV_LINE)
+		    n_zero_lines++;
+		else if (type & GV_BOUNDARY)
+		    n_zero_boundaries++;
+		    
+		if (Err)
+		    Vect_write_line(Err, type, Points, Cats);
+	    }
+	}
+    }
+    if (n_zero_lines)
+	G_warning(_("Number of lines of length zero: %d"), n_zero_lines);
+    if (n_zero_boundaries)
+	G_warning(_("Number of boundaries of length zero: %d"), n_zero_boundaries);
+
+    /* remaining checks are for areas only */
+    if (Vect_get_num_primitives(Map, GV_BOUNDARY) == 0)
+	return 1;
+
+    /* intersecting boundaries -> overlapping areas */
+    nerrors = Vect_check_line_breaks(Map, GV_BOUNDARY, Err);
+    if (nerrors)
+	G_warning(_("Number of boundary intersections: %d"), nerrors);
+
+    /* areas without centroids that are not holes
+     * only makes sense if all boundaries are correct */
+    nerrors = 0;
+    for (line = 1; line <= nlines; line++) {
+	int type;
+	
+	if (!Vect_line_alive(Map, line))
+	    continue;
+	    
+	type = Vect_get_line_type(Map, line);
+
+	if (type == GV_BOUNDARY) {
+	    struct P_topo_b *topo = (struct P_topo_b *)Map->plus.Line[line]->topo;
+
+	    if (topo->left == 0 || topo->right == 0) {
+		G_debug(3, "line = %d left = %d right = %d", line, 
+			topo->left, topo->right);
+		nerrors++;
+	    }
+	}
+    }
+    if (nerrors)
+	G_warning(_("Skipping further checks because of incorrect boundaries"));
+    else {
+	int i, area, left, right, neighbour;
+	int nareas = Vect_get_num_areas(Map);
+	struct ilist *List = Vect_new_list();
+
+	nerrors = 0;
+	for (area = 1; area <= nareas; area++) {
+	    if (!Vect_area_alive(Map, area))
+		continue;
+	    line = Vect_get_area_centroid(Map, area);
+	    if (line != 0)
+		continue;   /* has centroid */
+
+	    Vect_get_area_boundaries(Map, area, List);
+	    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;
+		    
+		if (neighbour < 0) {
+		    neighbour = Vect_get_isle_area(Map, abs(neighbour));
+		    if (!neighbour) {
+			/* borders outer void */
+			nerrors++;
+			if (Err) {
+			    Vect_read_line(Map, Points, Cats, abs(line));
+			    Vect_write_line(Err, GV_BOUNDARY, Points, Cats);
+			}
+		    }
+		    /* else neighbour is > 0, check below */
+		}
+		if (neighbour > 0) {
+		    if (Vect_get_area_centroid(Map, neighbour) == 0) {
+			/* neighbouring area does not have a centroid either */
+			nerrors++;
+			if (Err) {
+			    Vect_read_line(Map, Points, Cats, abs(line));
+			    Vect_write_line(Err, GV_BOUNDARY, Points, Cats);
+			}
+		    }
+		}
+	    }
+	}
+	Vect_destroy_list(List);
+
+	if (nerrors)
+	    G_warning(_("Number of redundant holes: %d"), 
+	              nerrors);
+    }
+
+    /* what else ? */
+
+    Vect_destroy_line_struct(Points);
+    Vect_destroy_cats_struct(Cats);
+
+    return 1;
+}
+
+/*!
    \brief Return current highest built level (part)
 
    \param Map vector map
@@ -642,19 +794,27 @@ int Vect_build_partial(struct Map_info *Map, int build)
 	G_message(_("Number of areas: %d"), plus->n_areas);
 	G_message(_("Number of isles: %d"), plus->n_isles);
 
+#if 0
+	/* not an error, message disabled to avoid confusion */
+	if (err_nocentr)
+	    G_message(_("Number of areas without centroid: %d"),
+		      err_nocentr);
+#endif
+
+	if (plus->n_clines > plus->n_areas)
+	    G_warning(_("Number of centroids exceeds number of areas: %d > %d"),
+		      plus->n_clines, plus->n_areas);
+
 	if (err_boundaries)
-	    G_message(_("Number of incorrect boundaries: %d"),
+	    G_warning(_("Number of incorrect boundaries: %d"),
 		      err_boundaries);
 
 	if (err_centr_out)
-	    G_message(_("Number of centroids outside area: %d"),
+	    G_warning(_("Number of centroids outside area: %d"),
 		      err_centr_out);
 
 	if (err_centr_dupl)
-	    G_message(_("Number of duplicate centroids: %d"), err_centr_dupl);
-
-	if (err_nocentr)
-	    G_message(_("Number of areas without centroid: %d"), err_nocentr);
+	    G_warning(_("Number of duplicate centroids: %d"), err_centr_dupl);
 
     }
     else if (build > GV_BUILD_NONE) {