Vlib: new Vect_snap_line()

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

include/defs/vector.h

@@ -354,8 +354,8 @@ int Vect_line_check_duplicate(const struct line_pnts *,
 void Vect_snap_lines(struct Map_info *, int, double, struct Map_info *);
 void Vect_snap_lines_list(struct Map_info *, const struct ilist *, double,
                           struct Map_info *);
-void Vect_snap_lines_list2(struct Map_info *, struct ilist *, int,
-                           double, struct Map_info *);
+int Vect_snap_line(struct Map_info *, struct ilist *, struct line_pnts *,
+                   double, int, int *, int *);
 void Vect_remove_dangles(struct Map_info *, int, double, struct Map_info *);
 void Vect_chtype_dangles(struct Map_info *, double, struct Map_info *);
 void Vect_select_dangles(struct Map_info *, int, double, struct ilist *);

lib/vector/Vlib/snap.c

@@ -26,12 +26,19 @@
 /* Vertex */
 typedef struct
 {
-    double x, y;
+    double x, y, z;
     int anchor;			/* 0 - anchor, do not snap this point, that means snap others to this */
     /* >0  - index of anchor to which snap this point */
     /* -1  - init value */
 } XPNT;
 
+/* Segment */
+typedef struct
+{
+    double x1, y1, z1,  /* start point */
+           x2, y2, z2;  /* end point */
+} XSEG;
+
 typedef struct
 {
     int anchor;
@@ -44,7 +51,7 @@ static int sort_new(const void *pa, const void *pb)
     NEW *p1 = (NEW *) pa;
     NEW *p2 = (NEW *) pb;
 
-    return (p1->along < p2->along ? -1 : 1);
+    return (p1->along < p2->along ? -1 : (p1->along > p2->along));
 
     /*
     if (p1->along < p2->along)
@@ -55,6 +62,20 @@ static int sort_new(const void *pa, const void *pb)
     */
 }
 
+typedef struct
+{
+    double x, y, z, along;
+} NEW2;
+
+/* for qsort */
+static int sort_new2(const void *pa, const void *pb)
+{
+    NEW2 *p1 = (NEW2 *) pa;
+    NEW2 *p2 = (NEW2 *) pb;
+
+    return (p1->along < p2->along ? -1 : (p1->along > p2->along));
+}
+
 /* This function is called by RTreeSearch() to find a vertex */
 static int find_item(int id, const struct RTree_Rect *rect, struct ilist *list)
 {
@@ -69,8 +90,42 @@ static int add_item(int id, const struct RTree_Rect *rect, struct ilist *list)
     return 1;
 }
 
+/* This function is called by RTreeSearch() to add selected node/line/area/isle to the list */
+static int find_item_box(int id, const struct RTree_Rect *rect, void *list)
+{
+    struct bound_box box;
+
+    box.W = rect->boundary[0];
+    box.S = rect->boundary[1];
+    box.B = rect->boundary[2];
+    box.E = rect->boundary[3];
+    box.N = rect->boundary[4];
+    box.T = rect->boundary[5];
+
+    dig_boxlist_add((struct boxlist *)list, id, &box);
+
+    return 0;
+}
+
+/* This function is called by RTreeSearch() to add selected node/line/area/isle to the list */
+static int add_item_box(int id, const struct RTree_Rect *rect, void *list)
+{
+    struct bound_box box;
+
+    box.W = rect->boundary[0];
+    box.S = rect->boundary[1];
+    box.B = rect->boundary[2];
+    box.E = rect->boundary[3];
+    box.N = rect->boundary[4];
+    box.T = rect->boundary[5];
+
+    dig_boxlist_add((struct boxlist *)list, id, &box);
+
+    return 1;
+}
+
 /* for ilist qsort'ing and bsearch'ing */
-int cmp_int(const void *a, const void *b)
+static int cmp_int(const void *a, const void *b)
 {
     int ai = *(int *)a;
     int bi = *(int *)b;
@@ -82,6 +137,7 @@ int cmp_int(const void *a, const void *b)
    \brief Snap selected lines to existing vertex in threshold.
    
    Snap selected lines to existing vertices of other selected lines.
+   3D snapping is not supported.
    
    Lines showing how vertices were snapped may be optionally written to error map. 
    Input map must be opened on level 2 for update at least on GV_BUILD_BASE.
@@ -475,462 +531,488 @@ Vect_snap_lines_list(struct Map_info *Map, const struct ilist *List_lines,
     G_verbose_message(_("New vertices: %d"), ncreated);
 }
 
+
 /*!
-   \brief Snap selected lines to existing vertex in threshold.
-   
-   Snap selected lines to existing vertices of lines of type type.
-   
-   Lines in the list will only be snapped to lines not in the list.
-   
-   If a list is given, new snapped lines are added to the list.
-   
-   If no list is given, all lines of type type are snapped to lines of 
-   type type.
-   
-   \param Map input map where vertices will be snapped
-   \param List_lines list of lines to snap
-   \param type type of lines where lines in list will be snapped to
-   \param thresh threshold in which snap vertices
+   \brief Snap lines in vector map to existing vertex in threshold.
+  
+   For details see Vect_snap_lines_list()
+  
+   \param[in] Map input map where vertices will be snapped
+   \param[in] type type of lines to snap
+   \param[in] thresh threshold in which snap vertices
    \param[out] Err vector map where lines representing snap are written or NULL
-   
+  
    \return void
-*/
+ */
 void
-Vect_snap_lines_list2(struct Map_info *Map, struct ilist *List_lines,
-		     int type, double thresh, struct Map_info *Err)
+Vect_snap_lines(struct Map_info *Map, int type, double thresh,
+		struct Map_info *Err)
 {
-    struct line_pnts *Points, *NPoints;
+    int line, nlines, ltype;
+    struct ilist *List;
+
+    List = Vect_new_list();
+
+    nlines = Vect_get_num_lines(Map);
+
+    for (line = 1; line <= nlines; line++) {
+	G_debug(3, "line =  %d", line);
+
+	if (!Vect_line_alive(Map, line))
+	    continue;
+
+	ltype = Vect_read_line(Map, NULL, NULL, line);
+
+	if (!(ltype & type))
+	    continue;
+
+	G_ilist_add(List, line);
+    }
+
+    Vect_snap_lines_list(Map, List, thresh, Err);
+
+    Vect_destroy_list(List);
+
+    return;
+}
+
+/*!
+   \brief Snap a line to reference lines in Map with threshold.
+  
+   3D snapping is supported. The line to snap and the reference lines 
+   can but do not need to be in different vector maps.
+   
+   Vect_snap_line() uses less memory, but is slower than 
+   Vect_snap_lines_list()
+
+   For details on snapping, see Vect_snap_lines_list()
+
+   \param[in] Map input map with reference lines
+   \param[in] reflist list of reference lines
+   \param[in,out] Points line points to snap
+   \param[in] thresh threshold in which to snap vertices
+   \param[in] with_z 2D or 3D snapping
+   \param[in,out] nsnapped number of snapped verices
+   \param[in,out] ncreated number of new vertices (on segments)
+  
+   \return 1 if line was changed, otherwise 0
+ */
+int
+Vect_snap_line(struct Map_info *Map, struct ilist *reflist,
+	       struct line_pnts *Points, double thresh, int with_z,
+	       int *nsnapped, int *ncreated)
+{
+    struct line_pnts *LPoints, *NPoints;
     struct line_cats *Cats;
-    int line, ltype, line_idx, nlines, new_line, val;
+    int i, v, line, nlines;
+    int changed;
     double thresh2;
 
     int point;			/* index in points array */
-    int nanchors, ntosnap;	/* number of anchors and number of points to be snapped */
-    int nsnapped, ncreated;	/* number of snapped verices, number of new vertices (on segments) */
-    int apoints, npoints, nvertices;	/* number of allocated points, registered points, vertices */
-    XPNT *XPnts;		/* Array of points */
-    NEW *New = NULL;		/* Array of new points */
+    int segment;		/* index in segments array */
+    int asegments;		/* number of allocated segments */
+    int nvertices;		/* number of vertices */
+    XSEG *XSegs = NULL;		/* Array of segments */
+    NEW2 *New = NULL;		/* Array of new points */
     int anew = 0, nnew;		/* allocated new points , number of new points */
-    struct ilist *List;
-    int *Index = NULL;		/* indexes of anchors for vertices */
-    int aindex = 0;		/* allocated Index */
+    struct boxlist *List;
 
-    struct RTree *RTree;
-    int rtreefd = -1;
+    struct RTree *pnt_tree,	/* spatial index for reference points */
+                 *seg_tree;	/* spatial index for reference segments */
     struct RTree_Rect rect;
 
     rect.boundary = G_malloc(6 * sizeof(RectReal));
-
-    if (List_lines) {
-	nlines = List_lines->n_values;
-	if (nlines < 1)
-	    return;
-	qsort(List_lines->value, nlines, sizeof(int), cmp_int);
-    }
-
-    Points = Vect_new_line_struct();
+    rect.boundary[0] = 0;
+    rect.boundary[1] = 0;
+    rect.boundary[2] = 0;
+    rect.boundary[3] = 0;
+    rect.boundary[4] = 0;
+    rect.boundary[5] = 0;
+
+    changed = 0;
+    if (nsnapped)
+	*nsnapped = 0;
+    if (ncreated)
+	*ncreated = 0;
+
+    point = Points->n_points;
+    Vect_line_prune(Points);
+    if (point != Points->n_points)
+	changed = 1;
+
+    nlines = reflist->n_values;
+    if (nlines < 1)
+	return changed;
+
+    LPoints = Vect_new_line_struct();
     NPoints = Vect_new_line_struct();
     Cats = Vect_new_cats_struct();
-    List = Vect_new_list();
-    if (getenv("GRASS_VECTOR_LOWMEM")) {
-	char *filename = G_tempfile();
-
-	rtreefd = open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
-	remove(filename);
-    }
-    RTree = RTreeCreateTree(rtreefd, 0, 2);
+    List = Vect_new_boxlist(1);
+    with_z = (with_z != 0);
+    pnt_tree = RTreeCreateTree(-1, 0, 2 + with_z);
+    RTreeSetOverflow(pnt_tree, 0);
+    seg_tree = RTreeCreateTree(-1, 0, 2 + with_z);
+    RTreeSetOverflow(seg_tree, 0);
 
     thresh2 = thresh * thresh;
 
-    /* Go through all lines in vector, and add each point to structure of points */
-    apoints = 0;
-    point = 1;			/* index starts from 1 ! */
+    point = segment = 1;	/* index starts from 1 ! */
     nvertices = 0;
-    XPnts = NULL;
+    asegments = 0;
 
-    nlines = Vect_get_num_lines(Map);
+    /* Add all vertices and all segments of all reference lines 
+     * to spatial indices */
+    nlines = reflist->n_values;
+    for (i = 0; i < nlines; i++) {
 
-    G_verbose_message(_("Snap vertices Pass 1: select points"));
-    for (line = 1; line <= nlines; line++) {
-	int v;
-
-	G_percent(line, nlines, 2);
+	line = reflist->value[i];
 
 	G_debug(3, "line =  %d", line);
 	if (!Vect_line_alive(Map, line))
 	    continue;
 
-	ltype = Vect_read_line(Map, Points, Cats, line);
-	val = -1;
-	if (List_lines) {
-	    if (bsearch(&line, List_lines->value, List_lines->n_values,
-			sizeof(int), cmp_int)) {
-		/* snap this line to other lines,
-		 * don't make its vertices anchor points */
-		val = -2;
-	    }
-	    else {
-		if (!(ltype & type))
-		    continue;
-	    }
-	}
-	else {
-	    if (!(ltype & type))
-		continue;
-	}
+	Vect_read_line(Map, LPoints, Cats, line);
+	Vect_line_prune(LPoints);
 
-	for (v = 0; v < Points->n_points; v++) {
+	for (v = 0; v < LPoints->n_points; v++) {
 	    G_debug(3, "  vertex v = %d", v);
 	    nvertices++;
 
 	    /* Box */
-	    rect.boundary[0] = Points->x[v];
-	    rect.boundary[3] = Points->x[v];
-	    rect.boundary[1] = Points->y[v];
-	    rect.boundary[4] = Points->y[v];
-	    rect.boundary[2] = 0;
-	    rect.boundary[5] = 0;
+	    rect.boundary[0] = LPoints->x[v];
+	    rect.boundary[3] = LPoints->x[v];
+	    rect.boundary[1] = LPoints->y[v];
+	    rect.boundary[4] = LPoints->y[v];
+	    if (with_z) {
+		rect.boundary[2] = LPoints->z[v];
+		rect.boundary[5] = LPoints->z[v];
+	    }
 
 	    /* Already registered ? */
-	    Vect_reset_list(List);
-	    RTreeSearch(RTree, &rect, (void *)find_item, List);
+	    Vect_reset_boxlist(List);
+	    RTreeSearch(pnt_tree, &rect, find_item_box, (void *)List);
 	    G_debug(3, "List : nvalues =  %d", List->n_values);
 
 	    if (List->n_values == 0) {	/* Not found */
-		/* Add to tree and to structure */
-		RTreeInsertRect(&rect, point, RTree);
-		if ((point - 1) == apoints) {
-		    apoints += 10000;
-		    XPnts =
-			(XPNT *) G_realloc(XPnts,
-					   (apoints + 1) * sizeof(XPNT));
-		}
-		XPnts[point].x = Points->x[v];
-		XPnts[point].y = Points->y[v];
-		XPnts[point].anchor = val;
-		point++;
-	    }
-	}
-    }
-
-    npoints = point - 1;
-
-    /* Go through all registered points and if not yet marked mark it as anchor and assign this anchor
-     * to all not yet marked points in threshold */
-
-    G_verbose_message(_("Snap vertices Pass 2: assign anchor vertices"));
-
-    nanchors = ntosnap = 0;
-    for (point = 1; point <= npoints; point++) {
-	int i;
-
-	G_percent(point, npoints, 2);
-
-	G_debug(3, "  point = %d", point);
-
-	/* also skip vertices of lines to snap */
-	if (XPnts[point].anchor >= 0 || XPnts[point].anchor == -2)
-	    continue;
 
-	XPnts[point].anchor = 0;	/* make it anchor */
-	nanchors++;
-
-	/* Find points in threshold */
-	rect.boundary[0] = XPnts[point].x - thresh;
-	rect.boundary[3] = XPnts[point].x + thresh;
-	rect.boundary[1] = XPnts[point].y - thresh;
-	rect.boundary[4] = XPnts[point].y + thresh;
-	rect.boundary[2] = 0;
-	rect.boundary[5] = 0;
-
-	Vect_reset_list(List);
-	RTreeSearch(RTree, &rect, (void *)add_item, List);
-	G_debug(4, "  %d points in threshold box", List->n_values);
+		/* Add to points tree */
+		RTreeInsertRect(&rect, point, pnt_tree);
 
-	for (i = 0; i < List->n_values; i++) {
-	    int pointb;
-	    double dx, dy, dist2;
-
-	    pointb = List->value[i];
-	    if (pointb == point)
-		continue;
-
-	    dx = XPnts[pointb].x - XPnts[point].x;
-	    dy = XPnts[pointb].y - XPnts[point].y;
-	    dist2 = dx * dx + dy * dy;
-
-	    if (dist2 > thresh2) /* outside threshold */
-		continue;
-		
-	    /* doesn't have an anchor yet */
-	    if (XPnts[pointb].anchor < 0) {
-		XPnts[pointb].anchor = point;
-		ntosnap++;
+		point++;
 	    }
-	    else if (XPnts[pointb].anchor > 0) {   /* check distance to previously assigned anchor */
-		double dist2_a;
+	    
+	    /* reference segments */
+	    if (v) {
+		/* Box */
+		if (LPoints->x[v - 1] < LPoints->x[v]) {
+		    rect.boundary[0] = LPoints->x[v - 1];
+		    rect.boundary[3] = LPoints->x[v];
+		}
+		else {
+		    rect.boundary[0] = LPoints->x[v];
+		    rect.boundary[3] = LPoints->x[v - 1];
+		}
+		if (LPoints->y[v - 1] < LPoints->y[v]) {
+		    rect.boundary[1] = LPoints->y[v - 1];
+		    rect.boundary[4] = LPoints->y[v];
+		}
+		else {
+		    rect.boundary[1] = LPoints->y[v];
+		    rect.boundary[4] = LPoints->y[v - 1];
+		}
+		if (LPoints->z[v - 1] < LPoints->z[v]) {
+		    rect.boundary[2] = LPoints->z[v - 1];
+		    rect.boundary[5] = LPoints->z[v];
+		}
+		else {
+		    rect.boundary[2] = LPoints->z[v];
+		    rect.boundary[5] = LPoints->z[v - 1];
+		}
 
-		dx = XPnts[XPnts[pointb].anchor].x - XPnts[pointb].x;
-		dy = XPnts[XPnts[pointb].anchor].y - XPnts[pointb].y;
-		dist2_a = dx * dx + dy * dy;
+		/* do not check for duplicates, too costly
+		 * because different segments can have identical boxes */
+		RTreeInsertRect(&rect, segment, seg_tree);
 
-		/* replace old anchor */
-		if (dist2 < dist2_a) {
-		    XPnts[pointb].anchor = point;
+		if ((segment - 1) == asegments) {
+		    asegments += 1000;
+		    XSegs =
+			(XSEG *) G_realloc(XSegs,
+					   (asegments + 1) * sizeof(XSEG));
+		}
+		XSegs[segment].x1 = LPoints->x[v - 1];
+		XSegs[segment].x2 = LPoints->x[v];
+		XSegs[segment].y1 = LPoints->y[v - 1];
+		XSegs[segment].y2 = LPoints->y[v];
+		if (with_z) {
+		    XSegs[segment].z1 = LPoints->z[v - 1];
+		    XSegs[segment].z1 = LPoints->z[v];
+		}
+		else {
+		    XSegs[segment].z1 = 0;
+		    XSegs[segment].z2 = 0;
 		}
+		segment++;
 	    }
 	}
     }
 
-    /* Go through all lines and: 
-     *   1) for all vertices: if not anchor snap it to its anchor
-     *   2) for all segments: snap it to all anchors in threshold (except anchors of vertices of course) */
-
-    nsnapped = ncreated = 0;
-
-    G_verbose_message(_("Snap vertices Pass 3: snap to assigned points"));
-
-    if (List_lines)
-	nlines = List_lines->n_values;
-    else
-	nlines = Vect_get_num_lines(Map);
-    for (line_idx = 0; line_idx < nlines; line_idx++) {
-	int v, spoint, anchor;
-	int changed = 0;
-
-	G_percent(line_idx, nlines, 2);
-
-	if (List_lines)
-	    line = List_lines->value[line_idx];
-	else
-	    line = line_idx;
-
-	G_debug(3, "line =  %d", line);
-	if (!Vect_line_alive(Map, line))
-	    continue;
-
-	ltype = Vect_read_line(Map, Points, Cats, line);
-
-	if (!List_lines && !(type & ltype))
-	    continue;
-
-	if (Points->n_points >= aindex) {
-	    aindex = Points->n_points;
-	    Index = (int *)G_realloc(Index, aindex * sizeof(int));
+    /* go through all vertices of the line to snap */
+    for (v = 0; v < Points->n_points; v++) {
+	double dist2, tmpdist2;
+	double x, y, z;
+
+	dist2 = thresh2 + thresh2;
+	x = Points->x[v];
+	y = Points->y[v];
+	z = Points->z[v];
+
+	/* Box */
+	rect.boundary[0] = Points->x[v] - thresh;
+	rect.boundary[3] = Points->x[v] + thresh;
+	rect.boundary[1] = Points->y[v] - thresh;
+	rect.boundary[4] = Points->y[v] + thresh;
+	if (with_z) {
+	    rect.boundary[2] = Points->z[v] - thresh;
+	    rect.boundary[5] = Points->z[v] + thresh;
 	}
 
-	/* Snap all vertices */
-	for (v = 0; v < Points->n_points; v++) {
-	    /* Box */
-	    rect.boundary[0] = Points->x[v];
-	    rect.boundary[3] = Points->x[v];
-	    rect.boundary[1] = Points->y[v];
-	    rect.boundary[4] = Points->y[v];
-	    rect.boundary[2] = 0;
-	    rect.boundary[5] = 0;
+	/* find nearest reference vertex */
+	Vect_reset_boxlist(List);
 
-	    /* Find point ( should always find one point ) */
-	    Vect_reset_list(List);
-
-	    RTreeSearch(RTree, &rect, (void *)add_item, List);
+	RTreeSearch(pnt_tree, &rect, add_item_box, (void *)List);
 
-	    spoint = List->value[0];
-	    anchor = XPnts[spoint].anchor;
-
-	    if (anchor > 0) {	/* to be snapped */
-		Points->x[v] = XPnts[anchor].x;
-		Points->y[v] = XPnts[anchor].y;
-		nsnapped++;
-		changed = 1;
-		Index[v] = anchor;	/* point on new location */
-	    }
-	    else {
-		Index[v] = spoint;	/* old point */
+	for (i = 0; i < List->n_values; i++) {
+	    double dx = List->box[i].E - Points->x[v];
+	    double dy = List->box[i].N - Points->y[v];
+	    double dz = 0;
+	    
+	    if (with_z)
+		dz = List->box[i].T - Points->z[v];
+	    
+	    tmpdist2 = dx * dx + dy * dy + dz * dz;
+	    
+	    if (tmpdist2 < dist2) {
+		dist2 = tmpdist2;
+
+		x = List->box[i].E;
+		y = List->box[i].N;
+		z = List->box[i].T;
 	    }
 	}
 
-	/* New points */
-	Vect_reset_line(NPoints);
+	if (dist2 <= thresh2 && dist2 > 0) {
+	    Points->x[v] = x;
+	    Points->y[v] = y;
+	    Points->z[v] = z;
 
-	/* Snap all segments to anchors in threshold */
-	for (v = 0; v < Points->n_points - 1; v++) {
-	    int i;
-	    double x1, x2, y1, y2, xmin, xmax, ymin, ymax;
+	    changed = 1;
+	    if (nsnapped)
+		(*nsnapped)++;
+	}
+    }
 
-	    G_debug(3, "  segment = %d end anchors : %d  %d", v, Index[v],
-		    Index[v + 1]);
+    /* go through all vertices of the line to snap */
+    for (v = 0; v < Points->n_points; v++) {
+	double dist2, tmpdist2;
+	double x, y, z;
+	
+	dist2 = thresh2 + thresh2;
+	x = Points->x[v];
+	y = Points->y[v];
+	z = Points->z[v];
+
+	/* Box */
+	rect.boundary[0] = Points->x[v] - thresh;
+	rect.boundary[3] = Points->x[v] + thresh;
+	rect.boundary[1] = Points->y[v] - thresh;
+	rect.boundary[4] = Points->y[v] + thresh;
+	if (with_z) {
+	    rect.boundary[2] = Points->z[v] - thresh;
+	    rect.boundary[5] = Points->z[v] + thresh;
+	}
 
-	    x1 = Points->x[v];
-	    x2 = Points->x[v + 1];
-	    y1 = Points->y[v];
-	    y2 = Points->y[v + 1];
+	/* find nearest reference segment */
+	Vect_reset_boxlist(List);
 
-	    Vect_append_point(NPoints, Points->x[v], Points->y[v],
-			      Points->z[v]);
+	RTreeSearch(seg_tree, &rect, add_item_box, (void *)List);
 
-	    /* Box */
-	    if (x1 <= x2) {
-		xmin = x1;
-		xmax = x2;
-	    }
-	    else {
-		xmin = x2;
-		xmax = x1;
-	    }
-	    if (y1 <= y2) {
-		ymin = y1;
-		ymax = y2;
-	    }
-	    else {
-		ymin = y2;
-		ymax = y1;
+	for (i = 0; i < List->n_values; i++) {
+	    double tmpx, tmpy, tmpz;
+	    int segment, status;
+	    
+	    segment = List->id[i];
+	    
+	    /* Check the distance */
+	    tmpdist2 =
+		dig_distance2_point_to_line(Points->x[v],
+					    Points->y[v],
+					    Points->z[v], 
+					    XSegs[segment].x1,
+					    XSegs[segment].y1,
+					    XSegs[segment].z1,
+					    XSegs[segment].x2,
+					    XSegs[segment].y2,
+					    XSegs[segment].z2,
+					    with_z, &tmpx, &tmpy, &tmpz,
+					    NULL, &status);
+
+	    if (tmpdist2 < dist2 && status == 0) {
+		dist2 = tmpdist2;
+
+		x = tmpx;
+		y = tmpy;
+		z = tmpz;
 	    }
+	}
 
-	    rect.boundary[0] = xmin - thresh;
-	    rect.boundary[3] = xmax + thresh;
-	    rect.boundary[1] = ymin - thresh;
-	    rect.boundary[4] = ymax + thresh;
-	    rect.boundary[2] = 0;
-	    rect.boundary[5] = 0;
-
-	    /* Find points */
-	    Vect_reset_list(List);
-	    RTreeSearch(RTree, &rect, (void *)add_item, List);
-
-	    G_debug(3, "  %d points in box", List->n_values);
+	if (dist2 <= thresh2 && dist2 > 0) {
+	    Points->x[v] = x;
+	    Points->y[v] = y;
+	    Points->z[v] = z;
+	    
+	    changed = 1;
+	    if (nsnapped)
+		(*nsnapped)++;
+	}
+    }
 
-	    /* Snap to anchor in threshold different from end points */
-	    nnew = 0;
-	    for (i = 0; i < List->n_values; i++) {
-		double dist2, along;
+    RTreeDestroyTree(seg_tree);
+    G_free(XSegs);
+
+    /* go through all segments of the line to snap */
+    /* find nearest reference vertex, add this vertex */
+    for (v = 0; v < Points->n_points - 1; v++) {
+	double x1, x2, y1, y2, z1, z2;
+	double xmin, xmax, ymin, ymax, zmin, zmax;
+
+	x1 = Points->x[v];
+	x2 = Points->x[v + 1];
+	y1 = Points->y[v];
+	y2 = Points->y[v + 1];
+	if (with_z) {
+	    z1 = Points->z[v];
+	    z2 = Points->z[v + 1];
+	}
+	else {
+	    z1 = z2 = 0;
+	}
 
-		spoint = List->value[i];
-		G_debug(4, "    spoint = %d anchor = %d", spoint,
-			XPnts[spoint].anchor);
+	Vect_append_point(NPoints, Points->x[v], Points->y[v],
+			  Points->z[v]);
 
-		if (spoint == Index[v] || spoint == Index[v + 1])
-		    continue;	/* end point */
-		if (XPnts[spoint].anchor > 0)
-		    continue;	/* point is not anchor */
+	/* Box */
+	if (x1 <= x2) {
+	    xmin = x1;
+	    xmax = x2;
+	}
+	else {
+	    xmin = x2;
+	    xmax = x1;
+	}
+	if (y1 <= y2) {
+	    ymin = y1;
+	    ymax = y2;
+	}
+	else {
+	    ymin = y2;
+	    ymax = y1;
+	}
+	if (z1 <= z2) {
+	    zmin = z1;
+	    zmax = z2;
+	}
+	else {
+	    zmin = z2;
+	    zmax = z1;
+	}
 
-		/* Check the distance */
-		dist2 =
-		    dig_distance2_point_to_line(XPnts[spoint].x,
-						XPnts[spoint].y, 0, x1, y1, 0,
-						x2, y2, 0, 0, NULL, NULL,
-						NULL, &along, NULL);
+	rect.boundary[0] = xmin - thresh;
+	rect.boundary[3] = xmax + thresh;
+	rect.boundary[1] = ymin - thresh;
+	rect.boundary[4] = ymax + thresh;
+	rect.boundary[2] = zmin - thresh;
+	rect.boundary[5] = zmax + thresh;
 
-		G_debug(4, "      distance = %lf", sqrt(dist2));
+	/* Find points */
+	Vect_reset_boxlist(List);
+	RTreeSearch(pnt_tree, &rect, add_item_box, (void *)List);
 
-		if (dist2 <= thresh2) {
-		    G_debug(4, "      anchor in thresh, along = %lf", along);
+	G_debug(3, "  %d points in box", List->n_values);
 
-		    if (nnew == anew) {
-			anew += 100;
-			New = (NEW *) G_realloc(New, anew * sizeof(NEW));
-		    }
-		    New[nnew].anchor = spoint;
-		    New[nnew].along = along;
-		    nnew++;
+	/* Snap to vertex in threshold different from end points */
+	nnew = 0;
+	for (i = 0; i < List->n_values; i++) {
+	    double dist2, along;
+	    int status;
+	    
+	    if (!with_z)
+		List->box[i].T = 0;
+
+	    if (Points->x[v] == List->box[i].E && 
+	        Points->y[v] == List->box[i].N &&
+		Points->z[v] == List->box[i].T)
+		continue;	/* start point */
+
+	    if (Points->x[v + 1] == List->box[i].E && 
+	        Points->y[v + 1] == List->box[i].N &&
+		Points->z[v + 1] == List->box[i].T)
+		continue;	/* end point */
+
+	    /* Check the distance */
+	    dist2 =
+		dig_distance2_point_to_line(List->box[i].E,
+					    List->box[i].N,
+					    List->box[i].T, x1, y1, z1,
+					    x2, y2, z2, with_z, NULL, NULL,
+					    NULL, &along, &status);
+
+	    if (dist2 <= thresh2 && status == 0) {
+		G_debug(4, "      anchor in thresh, along = %lf", along);
+
+		if (nnew == anew) {
+		    anew += 100;
+		    New = (NEW2 *) G_realloc(New, anew * sizeof(NEW2));
 		}
+		New[nnew].x = List->box[i].E;
+		New[nnew].y = List->box[i].N;
+		New[nnew].z = List->box[i].T;
+		New[nnew].along = along;
+		nnew++;
 	    }
-	    G_debug(3, "  nnew = %d", nnew);
-	    /* insert new vertices */
-	    if (nnew > 0) {
-		/* sort by distance along the segment */
-		qsort(New, sizeof(char) * nnew, sizeof(NEW), sort_new);
-
-		for (i = 0; i < nnew; i++) {
-		    anchor = New[i].anchor;
-		    /* Vect_line_insert_point ( Points, ++v, XPnts[anchor].x, XPnts[anchor].y, 0); */
-		    Vect_append_point(NPoints, XPnts[anchor].x,
-				      XPnts[anchor].y, 0);
-		    ncreated++;
-		}
-		changed = 1;
+	    G_debug(3, "dist: %g, thresh: %g", dist2, thresh2);
+	}
+	G_debug(3, "  nnew = %d", nnew);
+	/* insert new vertices */
+	if (nnew > 0) {
+	    /* sort by distance along the segment */
+	    qsort(New, nnew, sizeof(NEW2), sort_new2);
+
+	    for (i = 0; i < nnew; i++) {
+		Vect_append_point(NPoints, New[i].x, New[i].y, New[i].z);
+		if (ncreated)
+		    (*ncreated)++;
 	    }
+	    changed = 1;
 	}
+    }
 
-	/* append end point */
-	v = Points->n_points - 1;
-	Vect_append_point(NPoints, Points->x[v], Points->y[v], Points->z[v]);
+    /* append end point */
+    v = Points->n_points - 1;
+    Vect_append_point(NPoints, Points->x[v], Points->y[v], Points->z[v]);
 
-	if (changed) {		/* rewrite the line */
-	    Vect_line_prune(NPoints);	/* remove duplicates */
-	    if (NPoints->n_points > 1 || !(ltype & GV_LINES)) {
-		new_line = Vect_rewrite_line(Map, line, ltype, NPoints, Cats);
-		if (List_lines)
-		    G_ilist_add(List_lines, new_line);
-	    }
-	    else {
-		Vect_delete_line(Map, line);
-	    }
-	    if (Err) {
-		Vect_write_line(Err, ltype, Points, Cats);
-	    }
-	}
-    }				/* for each line */
-    G_percent(line_idx, nlines, 2); /* finish it */
+    if (Points->n_points != NPoints->n_points) {
+	Vect_line_prune(NPoints);	/* remove duplicates */
+	Vect_reset_line(Points);
+	Vect_append_points(Points, NPoints, GV_FORWARD);
+    }
 
-    Vect_destroy_line_struct(Points);
+    Vect_destroy_line_struct(LPoints);
     Vect_destroy_line_struct(NPoints);
     Vect_destroy_cats_struct(Cats);
-    G_free(XPnts);
-    G_free(Index);
+    Vect_destroy_boxlist(List);
     G_free(New);
-    RTreeDestroyTree(RTree);
-    if (rtreefd >= 0)
-	close(rtreefd);
-
+    RTreeDestroyTree(pnt_tree);
     G_free(rect.boundary);
 
-    G_verbose_message(_("Snapped vertices: %d"), nsnapped);
-    G_verbose_message(_("New vertices: %d"), ncreated);
-}
-
-
-/*!
-   \brief Snap lines in vector map to existing vertex in threshold.
-  
-   For details see Vect_snap_lines_list()
-  
-   \param[in] Map input map where vertices will be snapped
-   \param[in] type type of lines to snap
-   \param[in] thresh threshold in which snap vertices
-   \param[out] Err vector map where lines representing snap are written or NULL
-  
-   \return void
- */
-void
-Vect_snap_lines(struct Map_info *Map, int type, double thresh,
-		struct Map_info *Err)
-{
-    int line, nlines, ltype;
-    struct ilist *List;
-
-    List = Vect_new_list();
-
-    nlines = Vect_get_num_lines(Map);
-
-    for (line = 1; line <= nlines; line++) {
-	G_debug(3, "line =  %d", line);
-
-	if (!Vect_line_alive(Map, line))
-	    continue;
-
-	ltype = Vect_read_line(Map, NULL, NULL, line);
-
-	if (!(ltype & type))
-	    continue;
-
-	/* no need to check for duplicates:
-	 * use G_ilist_add() instead of Vect_list_append() */
-	G_ilist_add(List, line);
-    }
-
-    Vect_snap_lines_list(Map, List, thresh, Err);
-
-    Vect_destroy_list(List);
-
-    return;
+    return changed;
 }