v.distance: any feature to any feature

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

vector/v.distance/distance.c

@@ -0,0 +1,329 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include <grass/vector.h>
+#include "local_proto.h"
+
+/* TODO: geodesic distance for latlong */
+
+int get_line_box(const struct line_pnts *Points, 
+                 struct bound_box *box)
+{
+    int i;
+
+    if (Points->n_points == 0) {
+	box->E = box->W = box->N = box->S = box->T = box->B = 0.0 / 0.0;
+	return 0;
+    }
+
+    box->E = box->W = Points->x[0];
+    box->N = box->S = Points->y[0];
+    box->T = box->B = Points->z[0];
+
+    for (i = 1; i < Points->n_points; i++) {
+	if (box->E < Points->x[i])
+	    box->E = Points->x[i];
+	if (box->W > Points->x[i])
+	    box->W = Points->x[i];
+	if (box->N < Points->y[i])
+	    box->N = Points->y[i];
+	if (box->S > Points->y[i])
+	    box->S = Points->y[i];
+	if (box->T < Points->z[i])
+	    box->T = Points->z[i];
+	if (box->B > Points->z[i])
+	    box->B = Points->z[i];
+    }
+
+    return 1;
+}
+
+/* calculate distance parameters between two primitives
+ * return 1 point to point
+ * return 2 point to line
+ * return 1 line to line
+ */
+int line2line(struct line_pnts *FPoints, int ftype,
+              struct line_pnts *TPoints, int ttype,
+	      double *fx, double *fy, double *fz,
+	      double *falong, double *fangle,
+	      double *tx, double *ty, double *tz,
+	      double *talong, double *tangle,
+	      double *dist,
+	      int with_z,
+	      int geodesic)
+{
+    int i, fseg, tseg, tmp_seg;
+    double tmp_dist, tmp_x, tmp_y, tmp_z, tmp_along;
+    int ret = 1;
+    static struct line_pnts *iPoints = NULL;
+    
+    if (!iPoints)
+	iPoints = Vect_new_line_struct();
+
+    *dist = PORT_DOUBLE_MAX;
+
+    /* fangle and tangle are angles in radians, counter clockwise from x axis
+     * initialize to invalid angle */
+    *fangle = *tangle = -9.;
+    *falong = *talong = 0.;
+
+    *fx = FPoints->x[0];
+    *fy = FPoints->y[0];
+    *fz = FPoints->z[0];
+
+    *tx = TPoints->x[0];
+    *ty = TPoints->y[0];
+    *tz = TPoints->z[0];
+
+    /* point -> point */
+    if ((ftype & GV_POINTS) && (ttype & GV_POINTS)) {
+	Vect_line_distance(TPoints, FPoints->x[0], FPoints->y[0],
+			   FPoints->z[0], with_z, tx, ty,tz, dist, 
+			   NULL, talong);
+
+    }
+
+    /* point -> line and line -> line */
+    if ((ttype & GV_LINES)) {
+
+	tseg = 1;
+	/* calculate the min distance between each point in fline with tline */
+	for (i = 0; i < FPoints->n_points; i++) {
+
+	    tmp_seg = Vect_line_distance(TPoints, FPoints->x[i],
+	                                 FPoints->y[i], FPoints->z[i],
+					 with_z, &tmp_x, &tmp_y, &tmp_z,
+					 &tmp_dist, NULL, &tmp_along);
+	    if (*dist > tmp_dist) {
+		*dist = tmp_dist;
+		*tx = tmp_x;
+		*ty = tmp_y;
+		*tz = tmp_z;
+		*talong = tmp_along;
+		tseg = tmp_seg;
+	    }
+	}
+	Vect_point_on_line(TPoints, *talong, NULL, NULL, NULL,
+			   tangle, NULL);
+	ret++;
+    }
+
+    /* line -> point and line -> line */
+    if (ftype & GV_LINES) {
+	
+	fseg = 1;
+
+	/* calculate the min distance between each point in tline with fline */
+	for (i = 0; i < TPoints->n_points; i++) {
+
+	    tmp_seg = Vect_line_distance(FPoints, TPoints->x[i],
+			       TPoints->y[i], TPoints->z[i],
+			       with_z, &tmp_x, &tmp_y, &tmp_z,
+			       &tmp_dist, NULL, &tmp_along);
+	    if (*dist > tmp_dist) {
+		*dist = tmp_dist;
+		*fx = tmp_x;
+		*fy = tmp_y;
+		*fz = tmp_z;
+		*falong = tmp_along;
+		fseg = tmp_seg;
+	    }
+	}
+	Vect_point_on_line(FPoints, *falong, NULL, NULL, NULL,
+			   fangle, NULL);
+	ret++;
+	
+	if ((ttype & GV_LINES) && *dist > 0) {
+	    /* check for line intersection */
+	    struct bound_box fbox, tbox;
+
+	    get_line_box(FPoints, &fbox);
+	    get_line_box(TPoints, &tbox);
+
+	    if (Vect_box_overlap(&fbox, &tbox)) {
+		Vect_reset_line(iPoints);
+		Vect_line_get_intersections(FPoints, TPoints, iPoints, with_z);
+		if (iPoints->n_points) {
+		    *dist = 0;
+		    *fx = *tx = iPoints->x[0];
+		    *fy = *ty = iPoints->y[0];
+		    *fz = *tz = iPoints->z[0];
+		    
+		    /* falong, talong */
+		    Vect_line_distance(FPoints, iPoints->x[0],
+				       iPoints->y[0], iPoints->z[0],
+				       with_z, NULL, NULL, NULL,
+				       NULL, NULL, falong);
+		    Vect_line_distance(TPoints, iPoints->x[0],
+				       iPoints->y[0], iPoints->z[0],
+				       with_z, NULL, NULL, NULL,
+				       NULL, NULL, talong);
+		    /* fangle, tangle */
+		    Vect_point_on_line(FPoints, *falong, NULL, NULL, NULL,
+				       fangle, NULL);
+		    Vect_point_on_line(TPoints, *talong, NULL, NULL, NULL,
+				       tangle, NULL);
+		}
+	    }
+	}
+    }
+
+    return ret;
+}
+
+/* shortest distance between line and area
+ * return 1 inside area 
+ * return 2 inside isle of area 
+ * return 3 outside area */
+int line2area(const struct Map_info *To,
+	      struct line_pnts *Points, int type,
+	      int area, const struct bound_box *abox,
+	      double *fx, double *fy, double *fz,
+	      double *falong, double *fangle,
+	      double *tx, double *ty, double *tz,
+	      double *talong, double *tangle,
+	      double *dist,
+	      int with_z,
+	      int geodesic)
+{
+    int i, j;
+    double tmp_dist;
+    int isle, nisles;
+    static struct line_pnts *aPoints = NULL;
+    static struct line_pnts **iPoints = NULL;
+    static struct bound_box *ibox = NULL;
+    static int isle_alloc = 0;
+
+    if (!aPoints)
+	aPoints = Vect_new_line_struct();
+
+    *dist = PORT_DOUBLE_MAX;
+
+    /* fangle and tangle are angles in radians, counter clockwise from x axis
+     * initialize to invalid angle */
+    *fangle = *tangle = -9.;
+    *falong = *talong = 0.;
+
+    *fx = Points->x[0];
+    *fy = Points->y[0];
+    *fz = Points->z[0];
+
+    *tx = Points->x[0];
+    *ty = Points->y[0];
+    *tz = Points->z[0];
+
+    Vect_get_area_points(To, area, aPoints);
+    nisles = Vect_get_area_num_isles(To, area);
+    
+    if (nisles > isle_alloc) {
+	iPoints = G_realloc(iPoints, nisles * sizeof(struct line_pnts *));
+	ibox = G_realloc(ibox, nisles * sizeof(struct bound_box));
+	for (i = isle_alloc; i < nisles; i++)
+	    iPoints[i] = Vect_new_line_struct();
+	isle_alloc = nisles;
+    }
+    for (i = 0; i < nisles; i++) {
+	isle = Vect_get_area_isle(To, area, i);
+	Vect_get_isle_points(To, isle, iPoints[i]);
+	Vect_get_isle_box(To, isle, &ibox[i]);
+    }
+
+    /* inside area ? */
+
+    for (i = 0; i < Points->n_points; i++) {
+	if (Vect_point_in_box(Points->x[i], Points->y[i],
+			      Points->z[i], abox)) {
+
+	    int poly;
+	    
+	    poly = Vect_point_in_poly(Points->x[i], Points->y[i], aPoints);
+	    
+	    /* exactly on boundary */
+	    if (poly == 2) {
+		line2line(Points, type, aPoints, GV_BOUNDARY,
+		          fx, fy, fz, falong, fangle,
+		          tx, ty, tz, talong, tangle,
+			  dist, with_z, geodesic);
+		
+		return 1;
+	    }
+	    /* inside */
+	    else if (poly == 1) {
+		int inside = 0;
+
+		for (j = 0; j < nisles; j++) {
+		    if (Vect_point_in_box(Points->x[i], Points->y[i], Points->z[i], 
+					  &ibox[j])) {
+
+			poly = Vect_point_in_poly(Points->x[i], Points->y[i], iPoints[j]);
+			
+			/* inside or exactly on boundary */
+			if (poly > 0) {
+			    double tmp_fx, tmp_fy, tmp_fz, tmp_fangle, tmp_falong;
+			    double tmp_tx, tmp_ty, tmp_tz, tmp_tangle, tmp_talong;
+
+			    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,
+				      &tmp_dist, with_z, geodesic);
+
+			    if (*dist > tmp_dist) {
+				*dist = tmp_dist;
+
+				*fx = tmp_fx;
+				*fy = tmp_fy;
+				*fz = tmp_fz;
+				*falong = tmp_falong;
+				*fangle = tmp_fangle;
+
+				*tx = tmp_tx;
+				*ty = tmp_ty;
+				*tz = tmp_tz;
+				*talong = tmp_talong;
+				*tangle = tmp_tangle;
+			    }
+
+			    if (poly == 1)
+				inside = 1;
+			    
+			}
+		    }
+		    if (*dist == 0)
+			break;
+		}
+		/* inside area (inside outer ring, outside inner rings) */
+		if (!inside) {
+		    *fx = Points->x[i];
+		    *fy = Points->y[i];
+		    *fz = Points->z[i];
+
+		    *tx = Points->x[i];
+		    *ty = Points->y[i];
+		    *tz = Points->z[i];
+		    
+		    *dist = 0;
+
+		    return 1;
+		}
+		/* inside one of the islands */
+		else
+		    return 2;
+	    } /* end inside outer ring */
+	}
+	/* not possible */
+	if (*dist == 0)
+	    return 1;
+    }
+
+    /* line is outside area */
+
+    line2line(Points, type, aPoints, GV_BOUNDARY,
+	      fx, fy, fz, falong, fangle,
+	      tx, ty, tz, talong, tangle,
+	      dist, with_z, geodesic);
+
+    if (*dist == 0)
+	return 1;
+
+    return 3;
+}

vector/v.distance/local_proto.h

@@ -1,3 +1,4 @@
+#include <grass/vector.h>
 #include <grass/dbmi.h>
 
 /* define codes for characteristics of relation between two nearest features */
@@ -19,8 +20,10 @@ typedef struct
     int from_cat;		/* category (from) */
     int count;			/* number of features already found */
     int to_cat;			/* category (to) */
-    double from_x, from_y, from_z, to_x, to_y, to_z;	/* coordinates of nearest point */
+    double from_x, from_y, from_z; /* coordinates of nearest 'from' point */
+    double to_x, to_y, to_z;	/* coordinates of nearest 'to' point */
     double from_along, to_along;	/* distance along a linear feature to the nearest point */
+    double from_angle;		/* angle of linear feature in nearest point */
     double to_angle;		/* angle of linear feature in nearest point */
     double dist;		/* distance to nearest feature */
 } NEAR;
@@ -36,5 +39,29 @@ typedef struct
 int cmp_near(const void *, const void *);
 int cmp_near_to(const void *, const void *);
 int cmp_exist(const void *, const void *);
+
+/* distance.c */
+int get_line_box(const struct line_pnts *Points, 
+                 struct bound_box *box);
+int line2line(struct line_pnts *FPoints, int ftype,
+              struct line_pnts *TPoints, int ttype,
+	      double *fx, double *fy, double *fz,
+	      double *falong, double *fangle,
+	      double *tx, double *ty, double *tz,
+	      double *talong, double *tangle,
+	      double *dist,
+	      int with_z,
+	      int geodesic);
+int line2area(const struct Map_info *To,
+	      struct line_pnts *Points, int type,
+	      int area, const struct bound_box *abox,
+	      double *fx, double *fy, double *fz,
+	      double *falong, double *fangle,
+	      double *tx, double *ty, double *tz,
+	      double *talong, double *tangle,
+	      double *dist,
+	      int with_z,
+	      int geodesic);
+
 /* print.c */
 int print_upload(NEAR *, UPLOAD *, int, dbCatValArray *, dbCatVal *);

vector/v.distance/main.c

@@ -9,10 +9,11 @@
  *               - updated to 5.7 by Radim Blazek 2003
  *               - OGR support by Martin Landa <landa.martin gmail.com> (2009)
  *               - speed-up for dmax != 0 Markus Metz 2010
+ *               - support all features Markus Metz 2012
  *               
  * PURPOSE:      Calculates distance from a point to nearest feature in vector layer. 
  *               
- * COPYRIGHT:    (C) 2002-2010 by the GRASS Development Team
+ * COPYRIGHT:    (C) 2002-2012 by the GRASS Development Team
  *
  *               This program is free software under the 
  *               GNU General Public License (>=v2). 
@@ -34,7 +35,7 @@
 
 int main(int argc, char *argv[])
 {
-    int i, j, k;
+    int i, j;
     int print_as_matrix;	/* only for do_all=TRUE */
     int do_all;			/* calculate from each to each within the threshold */
     struct GModule *module;
@@ -49,7 +50,7 @@ int main(int argc, char *argv[])
     } flag;
     char *desc;
     struct Map_info From, To, Out, *Outp;
-    int from_type, to_type, from_field, to_field;
+    int from_type, to_type, from_field, to_field, with_z;
     double max, min;
     double *max_step;
     int n_max_steps, curr_step;
@@ -58,10 +59,12 @@ int main(int argc, char *argv[])
     NEAR *Near, *near;
     int anear;			/* allocated space, used only for do_all */
     UPLOAD *Upload;		/* zero terminated */
-    int ftype, fcat, tcat, count;
-    int nfrom, nto, nfcats, fline, tline, tseg, tarea, area, isle, nisles;
-    double tx, ty, tz, dist, talong, tmp_tx, tmp_ty, tmp_tz, tmp_dist,
-	tmp_talong;
+    int ftype, ttype, fcat, nfcats, tcat, count, fline, tfeature;
+    int nfrom, nfromlines, nfromareas, nto, ntolines, ntoareas;
+    int tarea, area, isle, nisles, nlines;
+    double fx, fy, fz, falong, fangle, tx, ty, tz, talong, tangle, dist;
+    double tmp_fx, tmp_fy, tmp_fz, tmp_falong, tmp_fangle, tmp_dist;
+    double tmp_tx, tmp_ty, tmp_tz, tmp_talong, tmp_tangle;
     struct field_info *Fi, *toFi;
     dbString stmt, dbstr;
     dbDriver *driver, *to_driver;
@@ -69,8 +72,8 @@ int main(int argc, char *argv[])
     char buf1[2000], buf2[2000], to_attr_sqltype[256];
     int update_ok, update_err, update_exist, update_notexist, update_dupl,
 	update_notfound, sqltype;
-    struct boxlist *List;
-    struct bound_box box;
+    struct boxlist *lList, *aList;
+    struct bound_box fbox, box;
     dbCatValArray cvarr;
     dbColumn *column;
 
@@ -100,8 +103,8 @@ int main(int argc, char *argv[])
 
     opt.from_type = G_define_standard_option(G_OPT_V_TYPE);
     opt.from_type->key = "from_type";
-    opt.from_type->options = "point,centroid";
-    opt.from_type->answer = "point";
+    opt.from_type->options = "point,line,boundary,centroid,area";
+    opt.from_type->answer = "point,line,area";
     opt.from_type->label = _("Feature type (from)");
     opt.from_type->guisection = _("From");
 
@@ -289,16 +292,36 @@ int main(int argc, char *argv[])
 
     from_field = Vect_get_field_number(&From, opt.from_field->answer);
 
-    if (Vect_get_num_primitives(&From, GV_POINTS) < 1) {
+    nfromlines = Vect_get_num_primitives(&From, from_type);
+    nfromareas = 0;
+    if (from_type & GV_AREA)
+	nfromareas = Vect_get_num_areas(&From);
+
+    nfrom = nfromlines + nfromareas;
+    if (nfrom < 1) {
 	const char *name = Vect_get_full_name(&From);
 	Vect_close(&From);
-	G_fatal_error(_("No points/centroids found in <%s>"), name);
+	G_fatal_error(_("No features of selected type found in <%s>"), name);
     }
     
     /* Open 'to' vector */
     Vect_set_open_level(2);
     Vect_open_old2(&To, opt.to->answer, "", opt.to_field->answer);
 
+    ntolines = Vect_get_num_primitives(&To, from_type);
+    ntoareas = 0;
+    if (to_type & GV_AREA)
+	ntoareas = Vect_get_num_areas(&To);
+
+    nto = ntolines + ntoareas;
+    if (nto < 1) {
+	const char *name = Vect_get_full_name(&To);
+	Vect_close(&From);
+	Vect_close(&To);
+	G_fatal_error(_("No features of selected type found in <%s>"), name);
+    }
+    with_z = (Vect_is_3d(&From) && Vect_is_3d(&To));
+
     to_field = Vect_get_field_number(&To, opt.to_field->answer);
 
     /* Open output vector */
@@ -317,7 +340,6 @@ int main(int argc, char *argv[])
     if (max != 0) {
 	struct bound_box fbox, tbox;
 	double dx, dy, dz, tmp_max;
-	int n_features = 0;
 
 	Vect_get_map_box(&From, &fbox);
 	Vect_get_map_box(&To, &tbox);
@@ -339,35 +361,24 @@ int main(int argc, char *argv[])
 	/* with max > 0 but max <<< tmp_max, 2 steps are sufficient, first 0 then max
 	 * a reasonable number of steps also depends on the number of features in To
 	 * e.g. only one area in To, no need to step */
-	nto = Vect_get_num_lines(&To);
-	for (tline = 1; tline <= nto; tline++) {
-	    /* TODO: Vect_get_line_type() */
-	    n_features += ((to_type & To.plus.Line[tline]->type) != 0);
-	}
-	if (to_type & GV_AREA) {
-	    if (Vect_get_num_areas(&To) > n_features)
-		n_features = Vect_get_num_areas(&To);
-	}
-	if (n_features == 0)
-	    G_fatal_error(_("No features of selected type in To vector <%s>"),
-			    opt.to->answer);
-	n_max_steps = sqrt(n_features) * max / tmp_max;
+
+	n_max_steps = sqrt(nto) * max / tmp_max;
 	/* max 9 steps from testing */
 	if (n_max_steps > 9)
 	    n_max_steps = 9;
 	if (n_max_steps < 2)
 	    n_max_steps = 2;
-	if (n_max_steps > n_features)
-	    n_max_steps = n_features;
+	if (n_max_steps > nto)
+	    n_max_steps = nto;
 
 	G_debug(2, "max = %f", max);
 	G_debug(2, "maximum reasonable search distance = %f", tmp_max);
-	G_debug(2, "n_features = %d", n_features);
+	G_debug(2, "n 'to' features = %d", nto);
 	G_debug(2, "n_max_steps = %d", n_max_steps);
     }
 
     if (min > max)
-	G_fatal_error("dmin can not be larger than dmax");
+	G_fatal_error(_("dmin can not be larger than dmax"));
 
     if (n_max_steps > 1) {
 	/* set up steps to increase search box */
@@ -391,6 +402,7 @@ int main(int argc, char *argv[])
     db_init_string(&stmt);
     db_init_string(&dbstr);
     driver = NULL;
+    Fi = NULL;
     if (!flag.print->answer && !do_all) {
 
 	Fi = Vect_get_field(&From, from_field);
@@ -424,6 +436,7 @@ int main(int argc, char *argv[])
     }
 
     to_driver = NULL;
+    toFi = NULL;
     if (opt.to_column->answer) {
 	toFi = Vect_get_field(&To, to_field);
 	if (toFi == NULL)
@@ -496,24 +509,19 @@ int main(int argc, char *argv[])
     TPoints = Vect_new_line_struct();
     FCats = Vect_new_cats_struct();
     TCats = Vect_new_cats_struct();
-    List = Vect_new_boxlist(1);
+    lList = Vect_new_boxlist(1); /* line list */
+    aList = Vect_new_boxlist(1); /* area list */
 
     /* Allocate space ( may be more than needed (duplicate cats and elements without cats) ) */
-    nfrom = Vect_get_num_lines(&From);
-    nto = Vect_get_num_lines(&To);
-    if (do_all) {
-	/* Be careful with do_all, it can easily run out of memory */
-	anear = 2 * nfrom;
-	Near = (NEAR *) G_calloc(anear, sizeof(NEAR));
-    }
-    else {
-	Near = (NEAR *) G_calloc(nfrom, sizeof(NEAR));
-    }
+    /* Be careful with do_all, it can easily run out of memory */
+    anear = nfrom;
+    Near = (NEAR *) G_calloc(nfrom, sizeof(NEAR));
 
     /* Read all cats from 'from' */
+    nfcats = 0;
     if (!do_all) {
-	nfcats = 0;
-	for (i = 1; i <= nfrom; i++) {
+	nlines = Vect_get_num_lines(&From);
+	for (i = 1; i <= nlines; i++) {
 	    ftype = Vect_read_line(&From, NULL, FCats, i);
 
 	    /* This keeps also categories of areas for future (if area s in from_type) */
@@ -525,22 +533,28 @@ int main(int argc, char *argv[])
 	    if (fcat < 0)
 		continue;
 	    Near[nfcats].from_cat = fcat;
+	    Near[nfcats].dist = -1;
+	    Near[nfcats].count = 0;
 	    nfcats++;
 	}
 	G_debug(1, "%d cats loaded from vector (including duplicates)",
 		nfcats);
+	
+	if (nfcats == 0)
+	    G_fatal_error(_("No categories for 'from' for slected type and layer"));
+
 	/* Sort by cats and remove duplicates */
 	qsort((void *)Near, nfcats, sizeof(NEAR), cmp_near);
 
 	/* remove duplicates */
+	j = 1;
 	for (i = 1; i < nfcats; i++) {
-	    if (Near[i].from_cat == Near[i - 1].from_cat) {
-		for (j = i; j < nfcats - 1; j++) {
-		    Near[j].from_cat = Near[j + 1].from_cat;
-		}
-		nfcats--;
+	    if (Near[i].from_cat != Near[j - 1].from_cat) {
+		Near[j].from_cat = Near[i].from_cat;
+		j++;
 	    }
 	}
+	nfcats = j;
 
 	G_debug(1, "%d cats loaded from vector (unique)", nfcats);
     }
@@ -549,20 +563,15 @@ int main(int argc, char *argv[])
     /* Note: as from_type is restricted to GV_POINTS (for now) everything is simple */
 
     count = 0;			/* count of distances in 'do_all' mode */
-    /* Find nearest lines */
-    if (to_type & (GV_POINTS | GV_LINES)) {
-	struct line_pnts *LLPoints;
+    /* Find nearest features for 'from' lines */
+    if (nfromlines) {
+	G_message(_("Finding nearest features..."));
+	
+	near = NULL;
+	nlines = Vect_get_num_lines(&From);
 
-	if (G_projection() == PROJECTION_LL) {
-	    LLPoints = Vect_new_line_struct();
-	}
-	else {
-	    LLPoints = NULL;
-	}
-	G_message(_("Finding nearest feature..."));
-	for (fline = 1; fline <= nfrom; fline++) {
+	for (fline = 1; fline <= nlines; fline++) {
 	    int tmp_tcat;
-	    double tmp_tangle, tangle;
 	    double tmp_min = (min < 0 ? 0 : min);
 	    double box_edge = 0;
 	    int done = FALSE;
@@ -571,99 +580,162 @@ int main(int argc, char *argv[])
 
 	    G_debug(3, "fline = %d", fline);
 	    G_percent(fline, nfrom, 2);
-	    ftype = Vect_read_line(&From, FPoints, FCats, fline);
+	    ftype = Vect_get_line_type(&From, fline);
 	    if (!(ftype & from_type))
 		continue;
 
+	    Vect_read_line(&From, FPoints, FCats, fline);
 	    Vect_cat_get(FCats, from_field, &fcat);
 	    if (fcat < 0 && !do_all)
 		continue;
 
+	    get_line_box(FPoints, &fbox);
+
+	    if (!do_all) {
+		/* find near by 'from' cat */
+		near = (NEAR *) bsearch((void *)&fcat, Near, nfcats,
+				        sizeof(NEAR), cmp_near);
+	    }
+
+	    dist = PORT_DOUBLE_MAX; /* distance to nearest 'to' feature */
+
 	    while (!done) {
 		done = TRUE;
 
+		tfeature = 0;  /* id of nearest 'to' feature */
+
 		if (!do_all) {
 		    /* enlarge search box until we get a hit */
 		    /* the objective is to enlarge the search box
 		     * in the first iterations just a little bit
 		     * to keep the number of hits low */
-		    Vect_reset_boxlist(List);
 		    while (curr_step < n_max_steps) {
 			box_edge = max_step[curr_step];
 
 			if (box_edge < tmp_min)
 			    continue;
 			
-			box.E = FPoints->x[0] + box_edge;
-			box.W = FPoints->x[0] - box_edge;
-			box.N = FPoints->y[0] + box_edge;
-			box.S = FPoints->y[0] - box_edge;
+			box.E = fbox.E + box_edge;
+			box.W = fbox.W - box_edge;
+			box.N = fbox.N + box_edge;
+			box.S = fbox.S - box_edge;
 			box.T = PORT_DOUBLE_MAX;
 			box.B = -PORT_DOUBLE_MAX;
 
-			Vect_select_lines_by_box(&To, &box, to_type, List);
+			if (ntolines)
+			    Vect_select_lines_by_box(&To, &box, to_type, lList);
+			if (ntoareas)
+			    Vect_select_areas_by_box(&To, &box, aList);
 
 			curr_step++;
-			if (List->n_values > 0)
+			if (lList->n_values > 0 || aList->n_values > 0)
 			    break;
 		    }
 		}
 		else {
-		    box.E = FPoints->x[0] + max;
-		    box.W = FPoints->x[0] - max;
-		    box.N = FPoints->y[0] + max;
-		    box.S = FPoints->y[0] - max;
+		    box.E = fbox.E + max;
+		    box.W = fbox.W - max;
+		    box.N = fbox.N + max;
+		    box.S = fbox.S - max;
 		    box.T = PORT_DOUBLE_MAX;
 		    box.B = -PORT_DOUBLE_MAX;
 
-		    Vect_select_lines_by_box(&To, &box, to_type, List);
+		    if (ntolines)
+			Vect_select_lines_by_box(&To, &box, to_type, lList);
+		    if (ntoareas)
+			Vect_select_areas_by_box(&To, &box, aList);
 		}
 
-		G_debug(3, "  %d lines in box", List->n_values);
+		G_debug(3, "  %d lines in box", lList->n_values);
 
-		tline = 0;
-		dist = PORT_DOUBLE_MAX;
-		for (i = 0; i < List->n_values; i++) {
+		for (i = 0; i < lList->n_values; i++) {
 		    tmp_tcat = -1;
-		    Vect_read_line(&To, TPoints, TCats, List->id[i]);
+		    ttype = Vect_read_line(&To, TPoints, TCats, lList->id[i]);
 
-		    tseg =
-			Vect_line_distance(TPoints, FPoints->x[0], FPoints->y[0],
-					   FPoints->z[0], (Vect_is_3d(&From) &&
-							   Vect_is_3d(&To)) ?
-					   WITH_Z : WITHOUT_Z, &tmp_tx, &tmp_ty,
-					   &tmp_tz, &tmp_dist, NULL, &tmp_talong);
-
-		    Vect_point_on_line(TPoints, tmp_talong, NULL, NULL, NULL,
-				       &tmp_tangle, NULL);
+		    line2line(FPoints, ftype, TPoints, ttype,
+		              &tmp_fx, &tmp_fy, &tmp_fz, &tmp_falong, &tmp_fangle,
+		              &tmp_tx, &tmp_ty, &tmp_tz, &tmp_talong, &tmp_tangle,
+			      &tmp_dist, with_z, 0);
 
 		    if (tmp_dist > max || tmp_dist < min)
 			continue;	/* not in threshold */
 
 		    /* TODO: more cats of the same field */
 		    Vect_cat_get(TCats, to_field, &tmp_tcat);
-		    if (G_projection() == PROJECTION_LL) {
-			/* calculate distances in meters not degrees (only 2D) */
-			Vect_reset_line(LLPoints);
-			Vect_append_point(LLPoints, FPoints->x[0], FPoints->y[0],
-					  FPoints->z[0]);
-			Vect_append_point(LLPoints, tmp_tx, tmp_ty, tmp_tz);
-			tmp_dist = Vect_line_geodesic_length(LLPoints);
-			Vect_reset_line(LLPoints);
-			for (k = 0; k < tseg; k++)
-			    Vect_append_point(LLPoints, TPoints->x[k],
-					      TPoints->y[k], TPoints->z[k]);
-			Vect_append_point(LLPoints, tmp_tx, tmp_ty, tmp_tz);
-			tmp_talong = Vect_line_geodesic_length(LLPoints);
+
+		    G_debug(4, "  tmp_dist = %f tmp_tcat = %d", tmp_dist,
+			    tmp_tcat);
+
+		    if (do_all) {
+			if (anear <= count) {
+			    anear += 10 + nto / 10;
+			    Near = (NEAR *) G_realloc(Near, anear * sizeof(NEAR));
+			    if (!Near)
+				G_fatal_error(_("Out of memory!"));
+			}
+			near = &(Near[count]);
+
+			/* store info about relation */
+			near->from_cat = fcat;
+			near->to_cat = tmp_tcat;	/* -1 is OK */
+			near->dist = tmp_dist;
+			near->from_x = tmp_fx;
+			near->from_y = tmp_fy;
+			near->from_z = tmp_fz;
+			near->from_along = tmp_falong;	/* 0 for points */
+			near->from_angle = tmp_fangle;
+			near->to_x = tmp_tx;
+			near->to_y = tmp_ty;
+			near->to_z = tmp_tz;
+			near->to_along = tmp_talong;	/* 0 for points */
+			near->to_angle = tmp_tangle;
+			near->count++;
+			count++;
+		    }
+		    else {
+			if (tfeature == 0 || (tmp_dist < dist)) {
+			    tfeature = lList->id[i];
+			    tcat = tmp_tcat;
+			    dist = tmp_dist;
+			    fx = tmp_fx;
+			    fy = tmp_fy;
+			    fz = tmp_fz;
+			    falong = tmp_falong;
+			    fangle = tmp_fangle;
+			    tx = tmp_tx;
+			    ty = tmp_ty;
+			    tz = tmp_tz;
+			    talong = tmp_talong;
+			    tangle = tmp_tangle;
+			}
 		    }
+		}
+
+		G_debug(3, "  %d areas in box", aList->n_values);
+
+		for (i = 0; i < aList->n_values; i++) {
+		    tmp_tcat = -1;
+
+		    line2area(&To, FPoints, ftype, aList->id[i], &aList->box[i],
+		              &tmp_fx, &tmp_fy, &tmp_fz, &tmp_falong, &tmp_fangle,
+		              &tmp_tx, &tmp_ty, &tmp_tz, &tmp_talong, &tmp_tangle,
+			      &tmp_dist, with_z, 0);
+
+		    if (tmp_dist > max || tmp_dist < min)
+			continue;	/* not in threshold */
+
+		    /* TODO: more cats of the same field */
+		    Vect_cat_get(TCats, to_field, &tmp_tcat);
 
 		    G_debug(4, "  tmp_dist = %f tmp_tcat = %d", tmp_dist,
 			    tmp_tcat);
 
 		    if (do_all) {
 			if (anear <= count) {
-			    anear += 10 + nfrom / 10;
+			    anear += 10 + nto / 10;
 			    Near = (NEAR *) G_realloc(Near, anear * sizeof(NEAR));
+			    if (!Near)
+				G_fatal_error(_("Out of memory!"));
 			}
 			near = &(Near[count]);
 
@@ -671,9 +743,11 @@ int main(int argc, char *argv[])
 			near->from_cat = fcat;
 			near->to_cat = tmp_tcat;	/* -1 is OK */
 			near->dist = tmp_dist;
-			near->from_x = FPoints->x[0];
-			near->from_y = FPoints->y[0];
-			near->from_z = FPoints->z[0];
+			near->from_x = tmp_fx;
+			near->from_y = tmp_fy;
+			near->from_z = tmp_fz;
+			near->from_along = tmp_falong;	/* 0 for points */
+			near->from_angle = tmp_fangle;
 			near->to_x = tmp_tx;
 			near->to_y = tmp_ty;
 			near->to_z = tmp_tz;
@@ -683,10 +757,15 @@ int main(int argc, char *argv[])
 			count++;
 		    }
 		    else {
-			if (tline == 0 || (tmp_dist < dist)) {
-			    tline = List->id[i];
+			if (tfeature == 0 || (tmp_dist < dist)) {
+			    tfeature = aList->id[i];
 			    tcat = tmp_tcat;
 			    dist = tmp_dist;
+			    fx = tmp_fx;
+			    fy = tmp_fy;
+			    fz = tmp_fz;
+			    falong = tmp_falong;
+			    fangle = tmp_fangle;
 			    tx = tmp_tx;
 			    ty = tmp_ty;
 			    tz = tmp_tz;
@@ -700,31 +779,29 @@ int main(int argc, char *argv[])
 
 		if (!do_all && curr_step < n_max_steps) {
 		    /* enlarging the search box is possible */
-		    if (tline > 0 && dist > box_edge) {
+		    if (tfeature > 0 && dist > box_edge) {
 			/* line found but distance > search edge:
 			 * line bbox overlaps with search box, line itself is outside search box */
 			done = FALSE;
 		    }
-		    else if (tline == 0) {
+		    else if (tfeature == 0) {
 			/* no line within max dist, but search box can still be enlarged */
 			done = FALSE;
 		    }
 		}
-		if (done && !do_all && tline > 0) {
-		    /* find near by cat */
-		    near =
-			(NEAR *) bsearch((void *)&fcat, Near, nfcats,
-					 sizeof(NEAR), cmp_near);
+		if (done && !do_all && tfeature > 0) {
 
-		    G_debug(4, "  near.from_cat = %d near.count = %d",
+		    G_debug(4, "  near->from_cat = %d near->count = %d",
 			    near->from_cat, near->count);
 		    /* store info about relation */
 		    if (near->count == 0 || near->dist > dist) {
 			near->to_cat = tcat;	/* -1 is OK */
 			near->dist = dist;
-			near->from_x = FPoints->x[0];
-			near->from_y = FPoints->y[0];
-			near->from_z = FPoints->z[0];
+			near->from_x = fx;
+			near->from_y = fy;
+			near->from_z = fz;
+			near->from_along = falong;	/* 0 for points */
+			near->from_angle = fangle;
 			near->to_x = tx;
 			near->to_y = ty;
 			near->to_z = tz;
@@ -733,117 +810,223 @@ int main(int argc, char *argv[])
 		    }
 		    near->count++;
 		}
-	    } /* done */
-	} /* next feature */
-	if (LLPoints) {
-	    Vect_destroy_line_struct(LLPoints);
-	}
+	    } /* done searching 'to' */
+	} /* next from feature */
     }
 
-    /* Find nearest areas */
-    if (to_type & GV_AREA) {
+    /* Find nearest features for 'from' areas */
+    /* the code is pretty much identical to the one for lines */
+    if (nfromareas) {
 	
-	G_message(_("Finding nearest areas..."));
-	for (fline = 1; fline <= nfrom; fline++) {
+	near = NULL;
+
+	G_message(_("Finding nearest features for areas..."));
+	for (area = 1; area <= nfromareas; area++) {
+	    int tmp_tcat;
 	    double tmp_min = (min < 0 ? 0 : min);
 	    double box_edge = 0;
 	    int done = FALSE;
 	    
 	    curr_step = 0;
 
-	    G_debug(3, "fline = %d", fline);
-	    G_percent(fline, nfrom, 2);
-	    ftype = Vect_read_line(&From, FPoints, FCats, fline);
-	    if (!(ftype & from_type))
-		continue;
+	    G_debug(3, "farea = %d", area);
+	    G_percent(area, nfromareas, 2);
 
+	    Vect_get_area_cats(&From, area, FCats);
 	    Vect_cat_get(FCats, from_field, &fcat);
 	    if (fcat < 0 && !do_all)
 		continue;
 
+	    Vect_get_area_box(&From, area, &fbox);
+	    Vect_reset_line(FPoints);
+
+	    if (!do_all) {
+		/* find near by 'from' cat */
+		near = (NEAR *) bsearch((void *)&fcat, Near, nfcats,
+				        sizeof(NEAR), cmp_near);
+	    }
+
+	    dist = PORT_DOUBLE_MAX; /* distance to nearest 'to' feature */
+
 	    while (!done) {
 		done = TRUE;
 
+		tfeature = 0;  /* id of nearest 'to' feature */
+
 		if (!do_all) {
 		    /* enlarge search box until we get a hit */
 		    /* the objective is to enlarge the search box
 		     * in the first iterations just a little bit
 		     * to keep the number of hits low */
-		    Vect_reset_boxlist(List);
 		    while (curr_step < n_max_steps) {
 			box_edge = max_step[curr_step];
 
 			if (box_edge < tmp_min)
 			    continue;
 			
-			box.E = FPoints->x[0] + box_edge;
-			box.W = FPoints->x[0] - box_edge;
-			box.N = FPoints->y[0] + box_edge;
-			box.S = FPoints->y[0] - box_edge;
+			box.E = fbox.E + box_edge;
+			box.W = fbox.W - box_edge;
+			box.N = fbox.N + box_edge;
+			box.S = fbox.S - box_edge;
 			box.T = PORT_DOUBLE_MAX;
 			box.B = -PORT_DOUBLE_MAX;
 
-			Vect_select_areas_by_box(&To, &box, List);
+			if (ntolines)
+			    Vect_select_lines_by_box(&To, &box, to_type, lList);
+			if (ntoareas)
+			    Vect_select_areas_by_box(&To, &box, aList);
 
 			curr_step++;
-			if (List->n_values > 0)
+			if (lList->n_values > 0 || aList->n_values > 0)
 			    break;
 		    }
 		}
 		else {
-		    box.E = FPoints->x[0] + max;
-		    box.W = FPoints->x[0] - max;
-		    box.N = FPoints->y[0] + max;
-		    box.S = FPoints->y[0] - max;
+		    box.E = fbox.E + max;
+		    box.W = fbox.W - max;
+		    box.N = fbox.N + max;
+		    box.S = fbox.S - max;
 		    box.T = PORT_DOUBLE_MAX;
 		    box.B = -PORT_DOUBLE_MAX;
 
-		    Vect_select_areas_by_box(&To, &box, List);
+		    if (ntolines)
+			Vect_select_lines_by_box(&To, &box, to_type, lList);
+		    if (ntoareas)
+			Vect_select_areas_by_box(&To, &box, aList);
 		}
 
-		G_debug(4, "%d areas selected by box", List->n_values);
+		G_debug(3, "  %d lines in box", lList->n_values);
+
+		for (i = 0; i < lList->n_values; i++) {
+		    tmp_tcat = -1;
+		    ttype = Vect_read_line(&To, TPoints, TCats, lList->id[i]);
+
+		    /* area to line */
+		    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 (tmp_dist > max || tmp_dist < min)
+			continue;	/* not in threshold */
+
+		    /* TODO: more cats of the same field */
+		    Vect_cat_get(TCats, to_field, &tmp_tcat);
+
+		    G_debug(4, "  tmp_dist = %f tmp_tcat = %d", tmp_dist,
+			    tmp_tcat);
+
+		    if (do_all) {
+			if (anear <= count) {
+			    anear += 10 + nto / 10;
+			    Near = (NEAR *) G_realloc(Near, anear * sizeof(NEAR));
+			    if (!Near)
+				G_fatal_error(_("Out of memory!"));
+			}
+			near = &(Near[count]);
+
+			/* store info about relation */
+			near->from_cat = fcat;
+			near->to_cat = tmp_tcat;	/* -1 is OK */
+			near->dist = tmp_dist;
+			near->from_x = tmp_fx;
+			near->from_y = tmp_fy;
+			near->from_z = tmp_fz;
+			near->from_along = tmp_falong;	/* 0 for points */
+			near->from_angle = tmp_fangle;
+			near->to_x = tmp_tx;
+			near->to_y = tmp_ty;
+			near->to_z = tmp_tz;
+			near->to_along = tmp_talong;	/* 0 for points */
+			near->to_angle = tmp_tangle;
+			near->count++;
+			count++;
+		    }
+		    else {
+			if (tfeature == 0 || (tmp_dist < dist)) {
+			    tfeature = lList->id[i];
+			    tcat = tmp_tcat;
+			    dist = tmp_dist;
+			    fx = tmp_fx;
+			    fy = tmp_fy;
+			    fz = tmp_fz;
+			    falong = tmp_falong;
+			    fangle = tmp_fangle;
+			    tx = tmp_tx;
+			    ty = tmp_ty;
+			    tz = tmp_tz;
+			    talong = tmp_talong;
+			    tangle = tmp_tangle;
+			}
+		    }
+		}
+
+		G_debug(3, "  %d areas in box", aList->n_values);
 
 		/* For each area in box check the distance */
-		tarea = 0;
-		dist = PORT_DOUBLE_MAX;
-		for (i = 0; i < List->n_values; i++) {
-		    int tmp_tcat;
+		for (i = 0; i < aList->n_values; i++) {
+		    int tmp_tcat, poly;
 
-		    area = List->id[i];
-		    G_debug(4, "%d: area %d", i, area);
-		    Vect_get_area_points(&To, area, TPoints);
+		    tarea = aList->id[i];
+		    G_debug(4, "%d: area %d", i, tarea);
+		    Vect_get_area_points(&To, tarea, TPoints);
+		    
+		    ttype = GV_BOUNDARY;
 
 		    /* Find the distance to this area */
-		    if (Vect_point_in_area(FPoints->x[0], FPoints->y[0], &To, area, &List->box[i])) {	/* in area */
-			tmp_dist = 0;
-			tmp_tx = FPoints->x[0];
-			tmp_ty = FPoints->y[0];
-		    }
-		    else if (Vect_point_in_poly(FPoints->x[0], FPoints->y[0], TPoints) > 0) {	/* in isle */
-			nisles = Vect_get_area_num_isles(&To, area);
-			for (j = 0; j < nisles; j++) {
-			    double tmp2_dist, tmp2_tx, tmp2_ty;
-
-			    isle = Vect_get_area_isle(&To, area, j);
-			    Vect_get_isle_points(&To, isle, TPoints);
-			    Vect_line_distance(TPoints, FPoints->x[0],
-					       FPoints->y[0], FPoints->z[0],
-					       WITHOUT_Z, &tmp2_tx, &tmp2_ty,
-					       NULL, &tmp2_dist, NULL, NULL);
-
-			    if (j == 0 || tmp2_dist < tmp_dist) {
-				tmp_dist = tmp2_dist;
-				tmp_tx = tmp2_tx;
-				tmp_ty = tmp2_ty;
+		    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,
+			 * check if 'from' area is inside 'to' area */
+			poly = 0;
+			/* boxes must overlap */
+			if (Vect_box_overlap(&fbox, &aList->box[i])) {
+			    if (FPoints->n_points == 0)
+				Vect_get_area_points(&From, area, FPoints);
+			    for (j = 0; j < FPoints->n_points; j++) {
+				poly = Vect_point_in_poly(FPoints->x[0], FPoints->y[0], TPoints);
+				if (poly)
+				    break;
+			    }
+			}
+			if (poly) {
+			    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;
+
+			    /* 'from' area is inside 'to' area,
+			     * get distance to 'to' isles */
+			    nisles = Vect_get_area_num_isles(&To, tarea);
+			    for (j = 0; j < nisles; j++) {
+				isle = Vect_get_area_isle(&To, tarea, j);
+				Vect_get_isle_points(&To, isle, TPoints);
+
+				line2area(&From, TPoints, ttype, area, &fbox,
+					  &tmp2_tx, &tmp2_ty, &tmp2_tz, &tmp2_talong, &tmp2_tangle,
+					  &tmp2_fx, &tmp2_fy, &tmp2_fz, &tmp2_falong, &tmp2_fangle,
+					  &tmp2_dist, with_z, 0);
+
+				if (tmp2_dist < tmp_dist) {
+				    tmp_dist = tmp2_dist;
+				    tmp_fx = tmp2_fx;
+				    tmp_fy = tmp2_fy;
+				    tmp_fz = tmp2_fz;
+				    tmp_falong = tmp2_falong;
+				    tmp_fangle = tmp2_fangle;
+				    tmp_tx = tmp2_tx;
+				    tmp_ty = tmp2_ty;
+				    tmp_tz = tmp2_tz;
+				    tmp_talong = tmp2_talong;
+				    tmp_tangle = tmp2_tangle;
+				}
 			    }
 			}
 		    }
-		    else {		/* outside area */
-			Vect_line_distance(TPoints, FPoints->x[0], FPoints->y[0],
-					   FPoints->z[0], WITHOUT_Z, &tmp_tx,
-					   &tmp_ty, NULL, &tmp_dist, NULL, NULL);
 
-		    }
 		    if (tmp_dist > max || tmp_dist < min)
 			continue;	/* not in threshold */
 		    Vect_get_area_cats(&To, area, TCats);
@@ -872,41 +1055,51 @@ int main(int argc, char *argv[])
 			near->from_cat = fcat;
 			near->to_cat = tmp_tcat;	/* -1 is OK */
 			near->dist = tmp_dist;
-			near->from_x = FPoints->x[0];
-			near->from_y = FPoints->y[0];
+			near->from_x = tmp_fx;
+			near->from_y = tmp_fy;
+			near->from_z = tmp_fz;
+			near->from_along = tmp_falong;	/* 0 for points */
+			near->from_angle = tmp_fangle;
 			near->to_x = tmp_tx;
 			near->to_y = tmp_ty;
-			near->to_along = 0;	/* nonsense for areas */
-			near->to_angle = 0;	/* not supported for areas */
+			near->to_z = tmp_tz;
+			near->to_along = tmp_talong;	/* 0 for points */
+			near->to_angle = tmp_tangle;
 			near->count++;
 			count++;
 		    }
-		    else if (tarea == 0 || tmp_dist < dist) {
-			tarea = area;
-			tcat = tmp_tcat;
-			dist = tmp_dist;
-			tx = tmp_tx;
-			ty = tmp_ty;
+		    else {
+			if (tfeature == 0 || tmp_dist < dist) {
+			    tfeature = tarea;
+			    tcat = tmp_tcat;
+			    dist = tmp_dist;
+			    fx = tmp_fx;
+			    fy = tmp_fy;
+			    fz = tmp_fz;
+			    falong = tmp_falong;
+			    fangle = tmp_fangle;
+			    tx = tmp_tx;
+			    ty = tmp_ty;
+			    tz = tmp_tz;
+			    talong = tmp_talong;
+			    tangle = tmp_tangle;
+			}
 		    }
 		}
 
 		if (!do_all && curr_step < n_max_steps) {
 		    /* enlarging the search box is possible */
-		    if (tarea > 0 && dist > box_edge) {
+		    if (tfeature > 0 && dist > box_edge) {
 			/* area found but distance > search edge:
 			 * area bbox overlaps with search box, area itself is outside search box */
 			done = FALSE;
 		    }
-		    else if (tarea == 0) {
+		    else if (tfeature == 0) {
 			/* no area within max dist, but search box can still be enlarged */
 			done = FALSE;
 		    }
 		}
-		if (done && !do_all && tarea > 0) {
-		    /* find near by cat */
-		    near =
-			(NEAR *) bsearch((void *)&fcat, Near, nfcats,
-					 sizeof(NEAR), cmp_near);
+		if (done && !do_all && tfeature > 0) {
 
 		    G_debug(4, "near.from_cat = %d near.count = %d dist = %f",
 			    near->from_cat, near->count, near->dist);
@@ -915,12 +1108,16 @@ int main(int argc, char *argv[])
 		    if (near->count == 0 || near->dist > dist) {
 			near->to_cat = tcat;	/* -1 is OK */
 			near->dist = dist;
-			near->from_x = FPoints->x[0];
-			near->from_y = FPoints->y[0];
+			near->from_x = fx;
+			near->from_y = fy;
+			near->from_z = fz;
+			near->from_along = falong;
+			near->from_angle = fangle;
 			near->to_x = tx;
 			near->to_y = ty;
-			near->to_along = 0;	/* nonsense for areas */
-			near->to_angle = 0;	/* not supported for areas */
+			near->to_z = tz;
+			near->to_along = talong;
+			near->to_angle = tangle;
 		    }
 		    near->count++;
 		}
@@ -968,6 +1165,8 @@ int main(int argc, char *argv[])
 	}
     }
 
+    update_ok = update_err = update_exist = update_notexist = update_dupl =
+	update_notfound = ncatexist = 0;
 
     /* Update database / print to stdout / create output map */
     if (flag.print->answer) {	/* print header */
@@ -1028,8 +1227,6 @@ int main(int argc, char *argv[])
 	    db_select_int(driver, Fi->table, Fi->key, NULL, &catexist);
 	G_debug(1, "%d cats selected from the table", ncatexist);
     }
-    update_ok = update_err = update_exist = update_notexist = update_dupl =
-	update_notfound = 0;
 
     if (!do_all) {
 	count = nfcats;

vector/v.distance/v.distance.html

@@ -1,23 +1,25 @@
 <h2>DESCRIPTION</h2>
 
-<em>v.distance</em> finds the nearest element in vector map
-(<em>to</em>) for elements in vector map (<em>from</em>). Various
-information about the vectors' relationships (distance, category, etc.) may be uploaded
-to the attribute table attached to the first vector map, or printed to
-'stdout'.  A new vector map may be created where lines connecting
-nearest points on features are written. <em>dmin</em> and/or <em>dmax</em> can be used to limit the search radius.
+<em>v.distance</em> finds the nearest element in vector map (<em>to</em>) 
+for elements in vector map (<em>from</em>). Various information about 
+the vectors' relationships (distance, category, etc.) may be uploaded to 
+the attribute table attached to the first vector map, or printed to 
+'stdout'. A new vector map may be created where lines connecting 
+nearest points on features are written. <em>dmin</em> and/or 
+<em>dmax</em> can be used to limit the search radius.
 
 <h2>NOTES</h2>
 
-If a nearest feature does not have a category, the attribute column is updated
-to <em>null</em>.  This is true also for areas, which means for example,
-that if a point is in an island (area WITHOUT category), <em>v.distance</em> 
-does not search for the nearest area WITH category; the island is identified 
-as the nearest and category updated to null.
-<p>The upload <em>column</em>(s) must already exist. Create one with <em>v.db.addcolumn</em>.
-<p>In lat-long locations <em>v.distance</em> gives distances (<em>dist</em>
-and <em>to_along</em>) in meters not in degrees calculated as geodesic
-distances on a sphere.
+If a nearest feature does not have a category, the attribute column is 
+updated to <em>null</em>.
+<p>The upload <em>column</em>(s) must already exist. Create one with 
+<em>v.db.addcolumn</em>.
+
+<!-- needs Vect_line_geodesic_distance()
+<p>In lat-long locations <em>v.distance</em> gives distances 
+(<em>dist</em> and <em>to_along</em>) in meters not in degrees 
+calculated as geodesic distances on a sphere.
+-->
 
 <h2>EXAMPLES</h2>
 
@@ -37,7 +39,7 @@ v.distance from=pnt to=ln upload=cat column=linecat
 For each point from vector map <b>pnt</b>, find the <em>nearest area</em>
 from map <b>ar</b> within the given threshold and write the related
 area categories to column <b>areacat</b> in an attribute table attached
-to vector map <b>pnt</b> (in the case that a point falls into a polygon area,
+to vector map <b>pnt</b> (in the case that a point falls into an area,
 the distance is zero):
 
 <div class="code"><pre>
@@ -143,6 +145,7 @@ Janne Soimasuo 1994, University of Joensuu, Faculty of Forestry, Finland<br>
 Cmd line coordinates support: Markus Neteler, ITC-irst, Trento, Italy<br>
 Updated for 5.1: Radim Blazek, ITC-irst, Trento, Italy<br>
 Matrix-like output by Martin Landa, FBK-irst, Trento, Italy<br>
-Improved processing speed: Markus Metz
+Improved processing speed: Markus Metz<br>
+Distance from any feature to any feature Markus Metz
 
 <p><i>Last changed: $Date$</i>