v.select: re-organize code to select features from vector map A by features from other vector map B (fixes https://trac.osgeo.org/grass/ticket/3361)

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

vector/v.select/geos.c

@@ -10,31 +10,31 @@
 static int relate_geos(struct Map_info *, const GEOSGeometry *,
 		       int, int, const char *, int);
 
-int line_relate_geos(struct Map_info *BIn, const GEOSGeometry * AGeom,
-		     int bline, int operator, const char *relate)
+int line_relate_geos(struct Map_info *AIn, const GEOSGeometry *BGeom,
+		     int aline, int operator, const char *relate)
 {
-    return relate_geos(BIn, AGeom, bline, operator, relate, 0);
+    return relate_geos(AIn, BGeom, aline, operator, relate, 0);
 }
 
-int area_relate_geos(struct Map_info *BIn, const GEOSGeometry * AGeom,
-		     int barea, int operator, const char *relate)
+int area_relate_geos(struct Map_info *AIn, const GEOSGeometry *BGeom,
+		     int aarea, int operator, const char *relate)
 {
-    return relate_geos(BIn, AGeom, barea, operator, relate, 1);
+    return relate_geos(AIn, BGeom, aarea, operator, relate, 1);
 }
 
-int relate_geos(struct Map_info *BIn, const GEOSGeometry * AGeom,
-		int bfid, int operator, const char *relate, int area)
+int relate_geos(struct Map_info *AIn, const GEOSGeometry *BGeom,
+		int afid, int operator, const char *relate, int area)
 {
-    GEOSGeometry *BGeom = NULL;
+    GEOSGeometry *AGeom = NULL;
     int found;
 
     found = 0;
     if (area)
-	BGeom = Vect_read_area_geos(BIn, bfid);
+	AGeom = Vect_read_area_geos(AIn, afid);
     else
-	BGeom = Vect_read_line_geos(BIn, bfid, NULL);
+	AGeom = Vect_read_line_geos(AIn, afid, NULL);
 
-    if (!BGeom)
+    if (!AGeom)
 	return 0;
 
     /* 
@@ -115,8 +115,8 @@ int relate_geos(struct Map_info *BIn, const GEOSGeometry * AGeom,
 	break;
     }
 
-    if (BGeom)
-	GEOSGeom_destroy(BGeom);
+    if (AGeom)
+	GEOSGeom_destroy(AGeom);
 
     return found;
 }

vector/v.select/main.c

@@ -113,9 +113,7 @@ int main(int argc, char *argv[])
 
     /* Alloc space for input lines array */
     ALines = (int *)G_calloc(Vect_get_num_lines(&(In[0])) + 1, sizeof(int));
-    AAreas = NULL;
-    if (flag.reverse->answer)
-	AAreas = (int *)G_calloc(Vect_get_num_areas(&(In[0])) + 1, sizeof(int));
+    AAreas = (int *)G_calloc(Vect_get_num_areas(&(In[0])) + 1, sizeof(int));
 
     /* Read field info */
     IFi = Vect_get_field(&(In[0]), ifield[0]);
@@ -150,46 +148,55 @@ int main(int argc, char *argv[])
     finishGEOS();
 #endif
 
-    if (nfound != 0) {
-
-    native = Vect_maptype(&Out) == GV_FORMAT_NATIVE;
-
-    nfields = Vect_cidx_get_num_fields(&(In[0]));
-    cats = (int **)G_malloc(nfields * sizeof(int *));
-    ncats = (int *)G_malloc(nfields * sizeof(int));
-    fields = (int *)G_malloc(nfields * sizeof(int));
-
-    /* Write lines */
-    if (!flag.table->answer && !native) {
-	/* Copy attributes for OGR output */
-	Vect_copy_map_dblinks(&(In[0]), &Out, TRUE);
-    }
-    
-    write_lines(&(In[0]), IFi, ALines, AAreas,
-		&Out, flag.table->answer ? 1 : 0, flag.reverse->answer ? 1 : 0,
-		nfields, fields, ncats, cats);
-
-    /* Copy tables */
-    if (!flag.table->answer && native) {
-	copy_tabs(&(In[0]), &Out,
-		  nfields, fields, ncats, cats);
+    if (!flag.reverse->answer) {
+	G_free(AAreas);
+	AAreas = NULL;
     }
 
-    /* print info about skipped features & close input maps */
-    for (iopt = 0; iopt < 2; iopt++) {
-        if (nskipped[iopt] > 0) {
-            G_warning(_("%d features from <%s> without category skipped"),
-                      nskipped[iopt], Vect_get_full_name(&(In[iopt])));
-        }
-        Vect_close(&(In[iopt]));
-    }
 
-    Vect_build(&Out);
-    Vect_close(&Out);
+    if (nfound != 0) {
 
-    G_done_msg(_("%d features written to output."), Vect_get_num_lines(&Out));
-    } else {
-    G_done_msg(_("No features found !"));
+	native = Vect_maptype(&Out) == GV_FORMAT_NATIVE;
+
+	nfields = Vect_cidx_get_num_fields(&(In[0]));
+	cats = (int **)G_malloc(nfields * sizeof(int *));
+	ncats = (int *)G_malloc(nfields * sizeof(int));
+	fields = (int *)G_malloc(nfields * sizeof(int));
+
+	/* Write lines */
+	if (!flag.table->answer && !native) {
+	    /* Copy attributes for OGR output */
+	    Vect_copy_map_dblinks(&(In[0]), &Out, TRUE);
+	}
+	
+	write_lines(&(In[0]), IFi, ALines, AAreas,
+		    &Out, flag.table->answer ? 1 : 0, flag.reverse->answer ? 1 : 0,
+		    nfields, fields, ncats, cats);
+
+	/* Copy tables */
+	if (!flag.table->answer && native) {
+	    copy_tabs(&(In[0]), &Out,
+		      nfields, fields, ncats, cats);
+	}
+
+	/* print info about skipped features & close input maps */
+	for (iopt = 0; iopt < 2; iopt++) {
+	    if (nskipped[iopt] > 0) {
+		G_warning(_("%d features from <%s> without category skipped"),
+			  nskipped[iopt], Vect_get_full_name(&(In[iopt])));
+	    }
+	    Vect_close(&(In[iopt]));
+	}
+
+	Vect_build(&Out);
+	Vect_close(&Out);
+
+	G_done_msg(_("%d features written to output."), Vect_get_num_lines(&Out));
+    }
+    else {
+	Vect_close(&In[0]);
+	Vect_close(&In[1]);
+	G_done_msg(_("No features found !"));
     }
 
     exit(EXIT_SUCCESS);

vector/v.select/overlap.c

@@ -40,43 +40,17 @@ void add_aarea(struct Map_info *In, int aarea, int *ALines, int *AAreas)
 
 /* Returns 1 if line1 from Map1 overlaps area2 from Map2,
  *         0 otherwise */
-int line_overlap_area(struct line_pnts *LPoints, struct Map_info *AMap,
-		      int area)
+int line_overlap_area(struct line_pnts *LPoints, struct line_pnts *OPoints,
+		      struct line_pnts **IPoints, int nisles)
 {
-    int i, nisles, isle;
-    static struct line_pnts *APoints = NULL;
-    static struct line_pnts **IPoints = NULL;
-    static int isles_alloc = 0;
-
-    G_debug(4, "line_overlap_area area = %d", area);
-
-    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();
-    }
-
-    Vect_get_area_points(AMap, area, APoints);
-    nisles = Vect_get_area_num_isles(AMap, area);
+    int i, isle;
 
-    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]);
-    }
+    G_debug(4, "line_overlap_area area");
 
     /* Try if any of line vertices is within area */
     for (i = 0; i < LPoints->n_points; i++) {
 
-	if (Vect_point_in_poly(LPoints->x[i], LPoints->y[i], APoints)) {
+	if (Vect_point_in_poly(LPoints->x[i], LPoints->y[i], OPoints)) {
 	    int inside = 1;
 	    
 	    for (isle = 0; isle < nisles; isle++) {
@@ -99,13 +73,12 @@ int line_overlap_area(struct line_pnts *LPoints, struct Map_info *AMap,
     /* Try intersections of line with area/isles boundary */
     /* Outer boundary */
 
-    if (Vect_line_check_intersection2(LPoints, APoints, 0)) {
+    if (Vect_line_check_intersection2(LPoints, OPoints, 0)) {
 	G_debug(4, "  -> line intersects outer area boundary");
 	return 1;
     }
 
     for (i = 0; i < nisles; i++) {
-
 	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,8 @@ int select_lines(struct Map_info *, int, int,
 
 /* overlap.c */
 void add_aarea(struct Map_info *, int, int *, int *);
-int line_overlap_area(struct line_pnts *, struct Map_info *, int);
+int line_overlap_area(struct line_pnts *, struct line_pnts *,
+                      struct line_pnts **, int);
 
 /* write.c */
 void write_lines(struct Map_info *, struct field_info *, int *, int *,

vector/v.select/select.c

@@ -7,336 +7,375 @@
 #include "proto.h"
 
 int select_lines(struct Map_info *aIn, int atype, int afield,
-                  struct Map_info *bIn, int btype, int bfield,
-                  int cat_flag, int operator, const char *relate,
-                  int *ALines, int *AAreas, int* nskipped)
+                 struct Map_info *bIn, int btype, int bfield,
+                 int cat_flag, int operator, const char *relate,
+                 int *ALines, int *AAreas, int* nskipped)
 {
-    int i;
-    int nalines, aline, ltype;
+    int i, ai;
+    int nblines, bline, ltype;
     int nfound = 0;
     
     struct line_pnts *APoints, *BPoints;
-    struct ilist *BoundList, *LList;
-    struct boxlist *List, *TmpList;
+    struct line_pnts *OPoints, **IPoints;
+    int isle, nisles, isles_alloc;
+    struct ilist *BoundList;
+    struct boxlist *List;
 
 #ifdef HAVE_GEOS
     initGEOS(G_message, G_fatal_error);
-    GEOSGeometry *AGeom = NULL;
+    GEOSGeometry *BGeom = NULL;
 #else
-    void *AGeom = NULL;
+    void *BGeom = NULL;
 #endif
 
     nskipped[0] = nskipped[1] = 0;
     APoints = Vect_new_line_struct();
     BPoints = Vect_new_line_struct();
+    OPoints = 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();
+    nisles = 0;
+    
     List = Vect_new_boxlist(1);
-    TmpList = Vect_new_boxlist(1);
     BoundList = Vect_new_list();
-    LList = Vect_new_list();
 
-    nalines = Vect_get_num_lines(aIn);
+    nblines = Vect_get_num_lines(bIn);
     
-    /* Lines in A. Go through all lines and mark those that meets condition */
-    if (atype & (GV_POINTS | GV_LINES)) {
+    /* Lines in B */
+    if (btype & (GV_POINTS | GV_LINES)) {
 	G_message(_("Processing features..."));
 	
-	G_percent(0, nalines, 2);
-	for (aline = 1; aline <= nalines; aline++) {
-	    struct bound_box abox;
+	G_percent(0, nblines, 2);
+	for (bline = 1; bline <= nblines; bline++) {
+	    struct bound_box bbox;
 
-	    G_debug(3, "aline = %d", aline);
-	    G_percent(aline, nalines, 2);	/* must be before any continue */
+	    G_debug(3, "bline = %d", bline);
+	    G_percent(bline, nblines, 2);	/* must be before any continue */
 
 	    /* Check category */
-	    if (!cat_flag && Vect_get_line_cat(aIn, aline, afield) < 0) {
-		nskipped[0]++;
+	    if (!cat_flag && Vect_get_line_cat(bIn, bline, bfield) < 0) {
+		nskipped[1]++;
 		continue;
 	    }
 
 	    /* Read line and check type */
 	    if (operator != OP_OVERLAP) {
 #ifdef HAVE_GEOS
-		AGeom = Vect_read_line_geos(aIn, aline, &ltype);
+		BGeom = Vect_read_line_geos(bIn, bline, &ltype);
 #endif
-		if (!(ltype & (GV_POINT | GV_LINE)))
-		    continue;
-
-		if (!AGeom)
+		if (!BGeom)
 		    G_fatal_error(_("Unable to read line id %d from vector map <%s>"),
-				  aline, Vect_get_full_name(aIn));
+				  bline, Vect_get_full_name(bIn));
 	    }
 	    else {
-		ltype = Vect_read_line(aIn, APoints, NULL, aline);
+		ltype = Vect_read_line(bIn, BPoints, NULL, bline);
 	    }
 	    
-	    if (!(ltype & atype))
+	    if (!(ltype & btype))
 		continue;
 	    
-	    Vect_get_line_box(aIn, aline, &abox);
+	    Vect_get_line_box(bIn, bline, &bbox);
 
-	    /* Check if this line overlaps any feature in B */
-	    /* x Lines in B */
-	    if (btype & (GV_POINTS | GV_LINES)) {
-		int found = 0;
+	    /* Check if this line overlaps any feature in A */
+	    /* x Lines in A */
+	    if (atype & (GV_POINTS | GV_LINES)) {
 		
 		/* Lines */
-		Vect_select_lines_by_box(bIn, &abox, btype, List);
-		for (i = 0; i < List->n_values; i++) {
-		    int bline;
+		Vect_select_lines_by_box(aIn, &bbox, atype, List);
+		for (ai = 0; ai < List->n_values; ai++) {
+		    int aline;
 		    
-		    bline = List->id[i];
-		    G_debug(3, "  bline = %d", bline);
+		    aline = List->id[ai];
+		    G_debug(3, "  aline = %d", aline);
+
+		    if (ALines[aline] == 1)
+			continue;
 		    
 		    /* Check category */
 		    if (!cat_flag &&
-			Vect_get_line_cat(bIn, bline, bfield) < 0) {
-			nskipped[1]++;
+			Vect_get_line_cat(aIn, aline, afield) < 0) {
+			nskipped[0]++;
 			continue;
 		    }
 		    
 		    if (operator != OP_OVERLAP) {
 #ifdef HAVE_GEOS
-			if(line_relate_geos(bIn, AGeom,
-					    bline, operator, relate)) {
-
-			    found = 1;
-			    break;
+			if (line_relate_geos(aIn, BGeom, aline,
+			                     operator, relate)) {
+			    ALines[aline] = 1;
+			    nfound += 1;
 			}
 #endif
 		    }
 		    else {
-			Vect_read_line(bIn, BPoints, NULL, bline);
+			Vect_read_line(aIn, APoints, NULL, aline);
 
-			if (Vect_line_check_intersection2(APoints, BPoints, 0)) {
-			    found = 1;
-			    break;
+			if (Vect_line_check_intersection2(BPoints, APoints, 0)) {
+			    ALines[aline] = 1;
+			    nfound += 1;
 			}
 		    }
 		}
-		
-		if (found) {
-		    ALines[aline] = 1;
-		    nfound += 1;
-		    continue;	/* Go to next A line */
-		}
 	    }
 	    
-	    /* x Areas in B. */
-	    if (btype & GV_AREA) {
+	    /* x Areas in A. */
+	    if (atype & GV_AREA) {
 		
-		Vect_select_areas_by_box(bIn, &abox, List);
-		for (i = 0; i < List->n_values; i++) {
-		    int barea;
+		Vect_select_areas_by_box(aIn, &bbox, List);
+		for (ai = 0; ai < List->n_values; ai++) {
+		    int aarea;
 		    
-		    barea = List->id[i];
-		    G_debug(3, "  barea = %d", barea);
+		    aarea = List->id[ai];
+		    G_debug(3, "  aarea = %d", aarea);
 		    
-		    if (Vect_get_area_cat(bIn, barea, bfield) < 0) {
-			nskipped[1]++;
+		    if (AAreas[aarea] == 1)
+			continue;
+		    
+		    if (Vect_get_area_centroid(aIn, aarea) < 1)
+			continue;
+
+		    if (!cat_flag &&
+		        Vect_get_area_cat(aIn, aarea, afield) < 0) {
+			nskipped[0]++;
 			continue;
 		    }
 
 		    if (operator != OP_OVERLAP) {
 #ifdef HAVE_GEOS
-			if(area_relate_geos(bIn, AGeom,
-					    barea, operator, relate)) {
-			    ALines[aline] = 1;
+			if (area_relate_geos(aIn, BGeom, aarea, 
+					     operator, relate)) {
+			    add_aarea(aIn, aarea, ALines, AAreas);
 			    nfound += 1;
-			    break;
 			}
 #endif
 		    }
 		    else {
-			if (line_overlap_area(APoints, bIn, barea)) {
-			    ALines[aline] = 1;
+			Vect_get_area_points(aIn, aarea, OPoints);
+			nisles = Vect_get_area_num_isles(aIn, aarea);
+			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(aIn, aarea, i);
+			    Vect_get_isle_points(aIn, isle, IPoints[i]);
+			}
+
+			if (line_overlap_area(BPoints, OPoints, IPoints, nisles)) {
+			    add_aarea(aIn, aarea, ALines, AAreas);
 			    nfound += 1;
-			    break;
 			}
 		    }
 		}
 	    }
-	    if (operator != OP_OVERLAP) {
 #ifdef HAVE_GEOS
-		GEOSGeom_destroy(AGeom);
-#endif
-		AGeom = NULL;
+	    if (operator != OP_OVERLAP) {
+		GEOSGeom_destroy(BGeom);
+		BGeom = NULL;
 	    }
+#endif
 	}
     }
     
-    /* Areas in A. */
-    if (atype & GV_AREA) {
-	int aarea, naareas;
+    /* Areas in B. */
+    if (btype & GV_AREA) {
+	int barea, nbareas, bcentroid;
 
 	G_message(_("Processing areas..."));
 	
-	naareas = Vect_get_num_areas(aIn);
+	nbareas = Vect_get_num_areas(bIn);
 
-	G_percent(0, naareas, 2);
-	for (aarea = 1; aarea <= naareas; aarea++) {
-	    struct bound_box abox;
+	G_percent(0, nbareas, 2);
+	for (barea = 1; barea <= nbareas; barea++) {
+	    struct bound_box bbox;
 
-	    G_percent(aarea, naareas, 1);
+	    G_percent(barea, nbareas, 1);
+
+	    if ((bcentroid = Vect_get_area_centroid(bIn, barea)) < 1)
+		continue;
 
-	    if (Vect_get_area_cat(aIn, aarea, afield) < 0) {
-		nskipped[0]++;
+	    if (!cat_flag &&
+	        Vect_get_area_cat(bIn, barea, bfield) < 0) {
+		nskipped[1]++;
 		continue;
 	    }
-	
-	    Vect_get_area_box(aIn, aarea, &abox);
-	    abox.T = PORT_DOUBLE_MAX;
-	    abox.B = -PORT_DOUBLE_MAX;
+
+	    Vect_read_line(bIn, BPoints, NULL, bcentroid);
+
+	    Vect_get_area_box(bIn, barea, &bbox);
+	    bbox.T = PORT_DOUBLE_MAX;
+	    bbox.B = -PORT_DOUBLE_MAX;
 
 	    if (operator != OP_OVERLAP) {
 #ifdef HAVE_GEOS
-		AGeom = Vect_read_area_geos(aIn, aarea);
+		BGeom = Vect_read_area_geos(bIn, barea);
 #endif
-		if (!AGeom)
+		if (!BGeom)
 		    G_fatal_error(_("Unable to read area id %d from vector map <%s>"),
-				  aarea, Vect_get_full_name(aIn));
+				  barea, Vect_get_full_name(bIn));
+	    }
+	    else {
+		Vect_get_area_points(bIn, barea, OPoints);
+		nisles = Vect_get_area_num_isles(bIn, barea);
+		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(bIn, barea, i);
+		    Vect_get_isle_points(bIn, isle, IPoints[i]);
+		}
 	    }
 
-	    /* x Lines in B */
-	    if (btype & (GV_POINTS | GV_LINES)) {
-		Vect_select_lines_by_box(bIn, &abox, btype, List);
+	    /* x Lines in A */
+	    if (atype & (GV_POINTS | GV_LINES)) {
+		Vect_select_lines_by_box(aIn, &bbox, atype, List);
 
-		for (i = 0; i < List->n_values; i++) {
-		    int bline;
+		for (ai = 0; ai < List->n_values; ai++) {
+		    int aline;
 
-		    bline = List->id[i];
+		    aline = List->id[ai];
+		    if (ALines[aline] == 1)
+			continue;
 
 		    if (!cat_flag &&
-			Vect_get_line_cat(bIn, bline, bfield) < 0) {
-			nskipped[1]++;
+			Vect_get_line_cat(aIn, aline, afield) < 0) {
+			nskipped[0]++;
 			continue;
 		    }
 		    
 		    if (operator != OP_OVERLAP) {
 #ifdef HAVE_GEOS
-			if(line_relate_geos(bIn, AGeom,
-					    bline, operator, relate)) {
-			    add_aarea(aIn, aarea, ALines, AAreas);
+			if (line_relate_geos(aIn, BGeom, aline,
+					     operator, relate)) {
+			    ALines[aline] = 1;
 			    nfound += 1;
-			    break;
 			}
 #endif
 		    }
 		    else {
-			Vect_read_line(bIn, BPoints, NULL, bline);
+			Vect_read_line(aIn, APoints, NULL, aline);
 
-			if (line_overlap_area(BPoints, aIn, aarea)) {
-			    add_aarea(aIn, aarea, ALines, AAreas);
+			if (line_overlap_area(APoints, OPoints, IPoints, nisles)) {
+			    ALines[aline] = 1;
 			    nfound += 1;
-			    continue;
 			}
 		    }
 		}
 	    }
 
-	    /* x Areas in B */
-	    if (btype & GV_AREA) {
-		int naisles;
-		int found = 0;
-
-		/* List of areas B */
+	    /* x Areas in A */
+	    if (atype & GV_AREA) {
 
-		/* Make a list of features forming area A */
-		Vect_reset_list(LList);
+		/* List of areas A */
+		Vect_select_areas_by_box(aIn, &bbox, List);
 
-		Vect_get_area_boundaries(aIn, aarea, BoundList);
-		for (i = 0; i < BoundList->n_values; i++) {
-		    Vect_list_append(LList, abs(BoundList->value[i]));
-		}
-
-		naisles = Vect_get_area_num_isles(aIn, aarea);
+		for (ai = 0; ai < List->n_values; ai++) {
+		    int found = 0;
+		    int aarea, acentroid;
 
-		for (i = 0; i < naisles; i++) {
-		    int j, aisle;
+		    aarea = List->id[ai];
+		    G_debug(3, "  aarea = %d", aarea);
+		    
+		    if (AAreas[aarea] == 1)
+			continue;
 
-		    aisle = Vect_get_area_isle(aIn, aarea, i);
+		    if ((acentroid = Vect_get_area_centroid(aIn, aarea)) < 1)
+			continue;
 
-		    Vect_get_isle_boundaries(aIn, aisle, BoundList);
-		    for (j = 0; j < BoundList->n_values; j++) {
-			Vect_list_append(LList, BoundList->value[j]);
+		    if (!cat_flag &&
+		        Vect_get_area_cat(aIn, aarea, afield) < 0) {
+			nskipped[0]++;
+			continue;
 		    }
-		}
 
-		Vect_select_areas_by_box(bIn, &abox, TmpList);
-
-		for (i = 0; i < LList->n_values; i++) {
-		    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;
-
-			barea = TmpList->id[j];
-			G_debug(3, "  barea = %d", barea);
-
-			if (Vect_get_area_cat(bIn, barea, bfield) < 0) {
-			    nskipped[1]++;
-			    continue;
+		    if (operator != OP_OVERLAP) {
+#ifdef HAVE_GEOS
+			if (area_relate_geos(aIn, BGeom, aarea,
+					     operator, relate)) {
+			    found = 1;
+			}
+#endif
+		    }
+		    else {
+			/* A inside B ? */
+			Vect_read_line(aIn, APoints, NULL, acentroid);
+			if (line_overlap_area(APoints, OPoints, IPoints, nisles)) {
+			    found = 1;
 			}
 
-			/* Check if any centroid of area B is in area A.
-			 * This test is important in if area B is completely within area A */
-			if ((bcentroid = Vect_get_area_centroid(bIn, barea)) > 0)
-			    Vect_read_line(bIn, BPoints, NULL, bcentroid);
-			else {
-			    double x, y;
+			/* B inside A ? */
+			if (!found) {
+			    struct bound_box abox;
 
-			    Vect_get_point_in_area(bIn, barea, &x, &y);
-			    Vect_reset_line(BPoints);
-			    Vect_append_point(BPoints, x, y, 0.0);
-			}
+			    Vect_get_area_box(aIn, aarea, &abox);
+			    abox.T = PORT_DOUBLE_MAX;
+			    abox.B = -PORT_DOUBLE_MAX;
 
-			if (operator != OP_OVERLAP) {
-#ifdef HAVE_GEOS
-			    if(area_relate_geos(bIn, AGeom,
-						barea, operator, relate)) {
+			    if (Vect_point_in_area(BPoints->x[0], BPoints->y[0], aIn,
+						   aarea, &abox)) {
 				found = 1;
-				break;
 			    }
-#endif
 			}
-			else {
-			    if (Vect_point_in_area(BPoints->x[0], BPoints->y[0], aIn,
-			                           aarea, &abox)) {
-				found = 1;
-				break;
+
+			/* A overlaps B ? */
+			if (!found) {
+			    Vect_get_area_boundaries(aIn, aarea, BoundList);
+			    for (i = 0; i < BoundList->n_values; i++) {
+				Vect_read_line(aIn, APoints, NULL, abs(BoundList->value[i]));
+				if (line_overlap_area(APoints, OPoints, IPoints, nisles)) {
+				    found = 1;
+				    break;
+				}
 			    }
-			    
-			    /* Check intersectin of lines from List with area B */
-			    if (line_overlap_area(APoints, bIn, barea)) {
-				found = 1;
-				break;
+			}
+
+			if (!found) {
+			    int j, naisles;
+
+			    naisles = Vect_get_area_num_isles(aIn, aarea);
+			    for (j = 0; j < naisles; j++) {
+
+				isle = Vect_get_area_isle(aIn, aarea, j);
+
+				Vect_get_isle_boundaries(aIn, isle, BoundList);
+				for (i = 0; i < BoundList->n_values; i++) {
+				    Vect_read_line(aIn, APoints, NULL, abs(BoundList->value[i]));
+				    if (line_overlap_area(APoints, OPoints, IPoints, nisles)) {
+					found = 1;
+					break;
+				    }
+				}
+				if (found)
+				    break;
 			    }
 			}
 		    }
 		    if (found) {
 			add_aarea(aIn, aarea, ALines, AAreas);
-		        nfound += 1;
-			break;
+			nfound += 1;
 		    }
 		}
 	    }
-	    if (operator != OP_OVERLAP) {
 #ifdef HAVE_GEOS
-		GEOSGeom_destroy(AGeom);
-#endif
-		AGeom = NULL;
+	    if (operator != OP_OVERLAP) {
+		GEOSGeom_destroy(BGeom);
+		BGeom = NULL;
 	    }
+#endif
 	}
     }
 
     Vect_destroy_line_struct(APoints);
     Vect_destroy_line_struct(BPoints);
     Vect_destroy_list(BoundList);
-    Vect_destroy_list(LList);
     Vect_destroy_boxlist(List);
-    Vect_destroy_boxlist(TmpList);
 
     return nfound;
 }

vector/v.select/write.c

@@ -6,7 +6,7 @@ void write_lines(struct Map_info *In, struct field_info *IFi, int *ALines,
 		 struct Map_info *Out, int table_flag, int reverse_flag,
 		 int nfields, int *fields, int *ncats, int **cats)
 {
-  int i, f, j, aline, nalines;
+    int i, f, j, aline, nalines;
     int atype;
     
     struct line_pnts *APoints;
@@ -59,14 +59,17 @@ void write_lines(struct Map_info *In, struct field_info *IFi, int *ALines,
 
 	if (!table_flag && (IFi != NULL)) {
 	    for (i = 0; i < ACats->n_cats; i++) {
+		f = -1;
 		for (j = 0; j < nfields; j++) {	/* find field */
 		    if (fields[j] == ACats->field[i]) {
 			f = j;
 			break;
 		    }
 		}
-		cats[f][ncats[f]] = ACats->cat[i];
-		ncats[f]++;
+		if (f >= 0) {
+		    cats[f][ncats[f]] = ACats->cat[i];
+		    ncats[f]++;
+		}
 	    }
 	}
     }