v.distance update for areas to areas

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

vector/v.distance/distance.c

@@ -197,6 +197,7 @@ int line2area(const struct Map_info *To,
     int i, j;
     double tmp_dist;
     int isle, nisles;
+    int all_inside_outer, all_outside_outer, all_outside_inner;
     static struct line_pnts *aPoints = NULL;
     static struct line_pnts **iPoints = NULL;
     static struct bound_box *ibox = NULL;
@@ -237,6 +238,8 @@ int line2area(const struct Map_info *To,
     }
 
     /* inside area ? */
+    all_inside_outer = all_outside_outer = 1;
+    all_outside_inner = 1;
 
     for (i = 0; i < Points->n_points; i++) {
 	if (Vect_point_in_box(Points->x[i], Points->y[i],
@@ -246,6 +249,15 @@ int line2area(const struct Map_info *To,
 	    
 	    poly = Vect_point_in_poly(Points->x[i], Points->y[i], aPoints);
 	    
+	    if (poly > 0) {
+		/* inside outer ring */
+		all_outside_outer = 0;
+	    }
+	    else {
+		/* outside outer ring */
+		all_inside_outer = 0;
+	    }
+	    
 	    /* exactly on boundary */
 	    if (poly == 2) {
 		line2line(Points, type, aPoints, GV_BOUNDARY,
@@ -255,9 +267,9 @@ int line2area(const struct Map_info *To,
 		
 		return 1;
 	    }
-	    /* inside */
+	    /* inside outer ring */
 	    else if (poly == 1) {
-		int inside = 0;
+		int inside_isle = 0;
 
 		for (j = 0; j < nisles; j++) {
 		    if (Vect_point_in_box(Points->x[i], Points->y[i], Points->z[i], 
@@ -270,6 +282,8 @@ int line2area(const struct Map_info *To,
 			    double tmp_fx, tmp_fy, tmp_fz, tmp_fangle, tmp_falong;
 			    double tmp_tx, tmp_ty, tmp_tz, tmp_tangle, tmp_talong;
 
+			    /* pass all points of the line, 
+			     * this will catch an intersection */
 			    line2line(Points, type, iPoints[j], GV_BOUNDARY,
 				      &tmp_fx, &tmp_fy, &tmp_fz, &tmp_falong, &tmp_fangle,
 				      &tmp_tx, &tmp_ty, &tmp_tz, &tmp_talong, &tmp_tangle,
@@ -291,16 +305,17 @@ int line2area(const struct Map_info *To,
 				*tangle = tmp_tangle;
 			    }
 
-			    if (poly == 1)
-				inside = 1;
+			    if (poly == 1) /* excludes isle boundary */
+				inside_isle = 1;
 			    
 			}
 		    }
 		    if (*dist == 0)
 			break;
 		}
-		/* inside area (inside outer ring, outside inner rings) */
-		if (!inside) {
+		/* inside area (inside outer ring, outside inner rings
+		 * or exactly on one of the inner rings) */
+		if (!inside_isle) {
 		    *fx = Points->x[i];
 		    *fy = Points->y[i];
 		    *fz = Points->z[i];
@@ -313,17 +328,48 @@ int line2area(const struct Map_info *To,
 
 		    return 1;
 		}
-		/* inside one of the islands */
-		else
-		    return 2;
+		else {
+		    /* inside one of the islands */
+		    all_outside_inner = 0;
+		    if (*dist == 0) {
+			/* the line intersected with the isle boundary
+			 * -> line is partially inside the area */
+			return 1;
+		    }
+		    /* else continue with next point */
+		}
 	    } /* end inside outer ring */
 	}
-	/* not possible */
+	else {
+	    /* point not in box of outer ring */
+	    all_inside_outer = 0;
+	}
+	/* exactly on boundary */
 	if (*dist == 0)
 	    return 1;
     }
-
-    /* line is outside area */
+    
+    /* if all points are inside the outer ring and inside inner rings,
+     * there could still be an intersection with one of the inner rings */
+     if (all_inside_outer) {
+	 if (all_outside_inner) {
+	     /* at least one point is really inside the area!
+	      * that should have been detected above */
+	     G_fatal_error(_("At least one point is really inside the area!"));
+	 }
+	 /* else all points are inside one of the area isles 
+	  * and we already have the minimum distance */
+	  return 2;
+     }
+
+    /* if at least one point was found to be inside the outer ring, 
+     * but no point really inside the area,
+     * and at least one point outside,
+     * then there must be an intersection of the line with both 
+     * the outer ring and one of the isle boundaries */
+
+    /* if all line points are outside of the area,
+     * intersection is still possible */
 
     line2line(Points, type, aPoints, GV_BOUNDARY,
 	      fx, fy, fz, falong, fangle,

vector/v.distance/main.c

@@ -968,19 +968,21 @@ int main(int argc, char *argv[])
 		    int tmp_tcat, poly;
 
 		    tarea = aList->id[i];
-		    G_debug(4, "%d: area %d", i, tarea);
+		    G_debug(4, "%d: 'to' area id %d", i, tarea);
+
 		    Vect_get_area_points(&To, tarea, TPoints);
 		    
 		    ttype = GV_BOUNDARY;
 
-		    /* Find the distance to this area */
+		    /* Find the distance of the outer ring of 'to' area
+		     * to 'from' area */
 		    poly = line2area(&From, TPoints, ttype, area, &fbox,
 		              &tmp_tx, &tmp_ty, &tmp_tz, &tmp_talong, &tmp_tangle,
 		              &tmp_fx, &tmp_fy, &tmp_fz, &tmp_falong, &tmp_fangle,
 			      &tmp_dist, with_z, 0);
 
 		    if (poly == 3) {
-			/* 'to' area is outside 'from' area,
+			/* 'to' outer ring is outside 'from' area,
 			 * check if 'from' area is inside 'to' area */
 			poly = 0;
 			/* boxes must overlap */
@@ -996,6 +998,8 @@ int main(int argc, char *argv[])
 			if (poly) {
 			    /* 'from' area is (partially) inside 'to' area,
 			     * get distance to 'to' area */
+			    if (FPoints->n_points == 0)
+				Vect_get_area_points(&From, area, FPoints);
 			    poly = line2area(&To, FPoints, ttype, tarea, &aList->box[i],
 				      &tmp_fx, &tmp_fy, &tmp_fz, &tmp_falong, &tmp_fangle,
 				      &tmp_tx, &tmp_ty, &tmp_tz, &tmp_talong, &tmp_tangle,
@@ -1004,7 +1008,7 @@ int main(int argc, char *argv[])
 			    /* inside isle ? */
 			    poly = poly == 2;
 			}
-			if (poly) {
+			if (poly == 1) {
 			    double tmp2_tx, tmp2_ty, tmp2_tz, tmp2_talong, tmp2_tangle;
 			    double tmp2_fx, tmp2_fy, tmp2_fz, tmp2_falong, tmp2_fangle;
 			    double tmp2_dist;
@@ -1040,7 +1044,7 @@ int main(int argc, char *argv[])
 
 		    if (tmp_dist > max || tmp_dist < min)
 			continue;	/* not in threshold */
-		    Vect_get_area_cats(&To, area, TCats);
+		    Vect_get_area_cats(&To, tarea, TCats);
 		    tmp_tcat = -1;
 		    /* TODO: all cats of given field ? */
 		    for (j = 0; j < TCats->n_cats; j++) {

vector/v.distance/v.distance.html

@@ -8,6 +8,38 @@ the attribute table attached to the first vector map, or printed to
 nearest points on features are written. <em>dmin</em> and/or 
 <em>dmax</em> can be used to limit the search radius.
 
+<p>
+For lines to lines, say line A to line B, <em>v.distance</em> calculates 
+the shortest distance of each vertex in A with each segment (not vertex) 
+in B. The module then calculates the shortest distance of each vertex in 
+B to each segment in A. The overall shortest distance of A points to B 
+segments and B points to A segments is used. Additionally, 
+<em>v.distance</em> checks for intersections. In case of intersections,
+the first intersection found is used and the distance set to zero.
+
+<p>
+For lines to areas, the distance is set to zero if a line is (partially) 
+inside an area. The first point of the line that is inside the area is 
+used as common point. The distance is also set to zero if the line 
+intersects with the outer ring or any of the inner rings (isles), in 
+which case the fist intersection is used as common point.
+
+<p>
+For areas to areas, the module checks first for overlap or if one area
+is (partially) inside the other area. This is computationally quite
+intensive. If the outer rings of the two areas do not overlap, the
+distance is calculated as above for lines to lines, treating the outer
+rings as two lines. Again, the first point encountered falling into an
+area is used as common point, or the first intersection point.
+
+<p>
+For anything else than points to lines, there can be several common
+locations with zero distance, and the common location would then be the 
+result of an overlay consisting of several points, lines, or areas.
+<em>v.distance</em> selects in these cases a single point, and does not 
+create an overlay like <em>v.overlay</em>.
+
+
 <h2>NOTES</h2>
 
 If a nearest feature does not have a category, the attribute column is