v.select: speed up OP_OVERLAP; Vlib: add Bentley-Ottmann algorithm to find line intersections (trunk https://trac.osgeo.org/grass/changeset/62091, https://trac.osgeo.org/grass/changeset/62045, https://trac.osgeo.org/grass/changeset/62089, partial https://trac.osgeo.org/grass/changeset/63830 for intersect2.c)

git-svn-id: https://svn.osgeo.org/grass/grass/branches/releasebranch_7_0@65423 15284696-431f-4ddb-bdfa-cd5b030d7da7

include/defs/vector.h

@@ -435,9 +435,16 @@ int Vect_line_intersection(struct line_pnts *, struct line_pnts *,
                            struct bound_box *, struct bound_box *,
                            struct line_pnts ***, struct line_pnts ***, int *,
                            int *, int);
+int Vect_line_intersection2(struct line_pnts *, struct line_pnts *,
+                            struct bound_box *, struct bound_box *,
+                            struct line_pnts ***, struct line_pnts ***, int *,
+                            int *, int);
 int Vect_line_check_intersection(struct line_pnts *, struct line_pnts *, int);
+int Vect_line_check_intersection2(struct line_pnts *, struct line_pnts *, int);
 int Vect_line_get_intersections(struct line_pnts *, struct line_pnts *,
                                 struct line_pnts *, int);
+int Vect_line_get_intersections2(struct line_pnts *, struct line_pnts *,
+                                 struct line_pnts *, int);
 char *Vect_subst_var(const char *, const struct Map_info *);
 
 /* Custom spatial index */

vector/v.select/overlap.c

@@ -38,28 +38,55 @@ void add_aarea(struct Map_info *In, int aarea, int *ALines)
 
 /* Returns 1 if line1 from Map1 overlaps area2 from Map2,
  *         0 otherwise */
-int line_overlap_area(struct Map_info *LMap, int line, struct Map_info *AMap,
-		      int area, struct bound_box box)
+int line_overlap_area(struct line_pnts *LPoints, struct Map_info *AMap,
+		      int area)
 {
     int i, nisles, isle;
-    static struct line_pnts *LPoints = NULL;
     static struct line_pnts *APoints = NULL;
+    static struct line_pnts **IPoints = NULL;
+    static int isles_alloc = 0;
 
-    G_debug(4, "line_overlap_area line = %d area = %d", line, area);
+    G_debug(4, "line_overlap_area area = %d", area);
 
-    if (!LPoints) {
-	LPoints = Vect_new_line_struct();
+    if (!APoints) {
 	APoints = Vect_new_line_struct();
+	isles_alloc = 10;
+	IPoints = G_malloc(isles_alloc * sizeof(struct line_pnts *));
+	for (i = 0; i < isles_alloc; i++)
+	    IPoints[i] = Vect_new_line_struct();
     }
 
-    /* Read line coordinates */
-    Vect_read_line(LMap, LPoints, NULL, line);
+    Vect_get_area_points(AMap, area, APoints);
+    nisles = Vect_get_area_num_isles(AMap, area);
+
+    if (nisles >= isles_alloc) {
+	IPoints = G_realloc(IPoints, (nisles + 10) * sizeof(struct line_pnts *));
+	for (i = isles_alloc; i < nisles + 10; i++)
+	    IPoints[i] = Vect_new_line_struct();
+	isles_alloc = nisles + 10;
+    }
+
+    for (i = 0; i < nisles; i++) {
+	isle = Vect_get_area_isle(AMap, area, i);
+	Vect_get_isle_points(AMap, isle, IPoints[i]);
+    }
 
     /* Try if any of line vertices is within area */
     for (i = 0; i < LPoints->n_points; i++) {
-	if (Vect_point_in_area(LPoints->x[i], LPoints->y[i], AMap, area, &box)) {
-	    G_debug(4, "  -> line vertex inside area");
-	    return 1;
+
+	if (Vect_point_in_poly(LPoints->x[i], LPoints->y[i], APoints)) {
+	    int inside = 1;
+	    
+	    for (isle = 0; isle < nisles; isle++) {
+		if (Vect_point_in_poly(LPoints->x[i], LPoints->y[i], IPoints[isle])) {
+		    inside = 0;
+		    break;
+		}
+	    }
+	    if (inside) {
+		G_debug(4, "  -> line vertex inside area");
+		return 1;
+	    }
 	}
     }
 
@@ -69,20 +96,15 @@ int line_overlap_area(struct Map_info *LMap, int line, struct Map_info *AMap,
 
     /* Try intersections of line with area/isles boundary */
     /* Outer boundary */
-    Vect_get_area_points(AMap, area, APoints);
 
-    if (Vect_line_check_intersection(LPoints, APoints, 0)) {
+    if (Vect_line_check_intersection2(LPoints, APoints, 0)) {
 	G_debug(4, "  -> line intersects outer area boundary");
 	return 1;
     }
 
-    nisles = Vect_get_area_num_isles(AMap, area);
-
     for (i = 0; i < nisles; i++) {
-	isle = Vect_get_area_isle(AMap, area, i);
-	Vect_get_isle_points(AMap, isle, APoints);
 
-	if (Vect_line_check_intersection(LPoints, APoints, 0)) {
+	if (Vect_line_check_intersection2(LPoints, IPoints[i], 0)) {
 	    G_debug(4, "  -> line intersects area island boundary");
 	    return 1;
 	}

vector/v.select/proto.h

@@ -42,7 +42,7 @@ void select_lines(struct Map_info *, int, int,
 
 /* overlap.c */
 void add_aarea(struct Map_info *, int, int *);
-int line_overlap_area(struct Map_info *, int, struct Map_info *, int, struct bound_box);
+int line_overlap_area(struct line_pnts *, struct Map_info *, int);
 
 /* write.c */
 void write_lines(struct Map_info *, struct field_info *, int *,

vector/v.select/select.c

@@ -34,15 +34,12 @@ void select_lines(struct Map_info *aIn, int atype, int afield,
     LList = Vect_new_list();
 
     nalines = Vect_get_num_lines(aIn);
-    if (operator == OP_OVERLAP)
-	G_message("Using GRASS GIS operator...");
-    else
-	G_message("Using GEOS operator...");
     
     /* Lines in A. Go through all lines and mark those that meets condition */
     if (atype & (GV_POINTS | GV_LINES)) {
 	G_message(_("Processing features..."));
 	
+	G_percent(0, nalines, 2);
 	for (aline = 1; aline <= nalines; aline++) {
 	    struct bound_box abox;
 
@@ -109,7 +106,7 @@ void select_lines(struct Map_info *aIn, int atype, int afield,
 		    else {
 			Vect_read_line(bIn, BPoints, NULL, bline);
 
-			if (Vect_line_check_intersection(APoints, BPoints, 0)) {
+			if (Vect_line_check_intersection2(APoints, BPoints, 0)) {
 			    found = 1;
 			    break;
 			}
@@ -147,7 +144,7 @@ void select_lines(struct Map_info *aIn, int atype, int afield,
 #endif
 		    }
 		    else {
-			if (line_overlap_area(aIn, aline, bIn, barea, List->box[i])) {
+			if (line_overlap_area(APoints, bIn, barea)) {
 			    ALines[aline] = 1;
 			    break;
 			}
@@ -171,10 +168,11 @@ void select_lines(struct Map_info *aIn, int atype, int afield,
 	
 	naareas = Vect_get_num_areas(aIn);
 
+	G_percent(0, naareas, 2);
 	for (aarea = 1; aarea <= naareas; aarea++) {
 	    struct bound_box abox;
 
-	    G_percent(aarea, naareas, 2);	/* must be before any continue */
+	    G_percent(aarea, naareas, 1);
 
 	    if (Vect_get_area_cat(aIn, aarea, afield) < 0) {
 		nskipped[0]++;
@@ -219,7 +217,9 @@ void select_lines(struct Map_info *aIn, int atype, int afield,
 #endif
 		    }
 		    else {
-			if (line_overlap_area(bIn, bline, aIn, aarea, abox)) {
+			Vect_read_line(bIn, BPoints, NULL, bline);
+
+			if (line_overlap_area(BPoints, aIn, aarea)) {
 			    add_aarea(aIn, aarea, ALines);
 			    continue;
 			}
@@ -261,6 +261,7 @@ void select_lines(struct Map_info *aIn, int atype, int afield,
 		    int j;
 
 		    aline = abs(LList->value[i]);
+		    Vect_read_line(aIn, APoints, NULL, aline);
 
 		    for (j = 0; j < TmpList->n_values; j++) {
 			int barea, bcentroid;
@@ -302,8 +303,7 @@ void select_lines(struct Map_info *aIn, int atype, int afield,
 			    }
 			    
 			    /* Check intersectin of lines from List with area B */
-			    if (line_overlap_area(aIn, aline,
-						  bIn, barea, TmpList->box[j])) {
+			    if (line_overlap_area(APoints, bIn, barea)) {
 				found = 1;
 				break;
 			    }