[GRASS-SVN] r54135 - in grass/branches/releasebranch_6_4: include lib/vector/Vlib

svn_grass at osgeo.org svn_grass at osgeo.org
Sat Dec 1 12:30:22 PST 2012


Author: mmetz
Date: 2012-12-01 12:30:22 -0800 (Sat, 01 Dec 2012)
New Revision: 54135

Modified:
   grass/branches/releasebranch_6_4/include/Vect.h
   grass/branches/releasebranch_6_4/lib/vector/Vlib/snap.c
Log:
Vlib: new Vect_snap_line()

Modified: grass/branches/releasebranch_6_4/include/Vect.h
===================================================================
--- grass/branches/releasebranch_6_4/include/Vect.h	2012-12-01 20:03:35 UTC (rev 54134)
+++ grass/branches/releasebranch_6_4/include/Vect.h	2012-12-01 20:30:22 UTC (rev 54135)
@@ -331,6 +331,8 @@
 void Vect_snap_lines(struct Map_info *, int, double, struct Map_info *);
 void Vect_snap_lines_list(struct Map_info *, struct ilist *, double,
 			  struct Map_info *);
+int Vect_snap_line(struct Map_info *, struct ilist *, struct line_pnts *,
+                   double, 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 *);

Modified: grass/branches/releasebranch_6_4/lib/vector/Vlib/snap.c
===================================================================
--- grass/branches/releasebranch_6_4/lib/vector/Vlib/snap.c	2012-12-01 20:03:35 UTC (rev 54134)
+++ grass/branches/releasebranch_6_4/lib/vector/Vlib/snap.c	2012-12-01 20:30:22 UTC (rev 54135)
@@ -36,19 +36,47 @@
     /* -1  - init value */
 } XPNT;
 
+/* Segment */
 typedef struct
 {
+    double x1, y1,  /* start point */
+           x2, y2;  /* end point */
+} XSEG;
+
+typedef struct
+{
     int anchor;
     double along;
 } NEW;
 
-/* This function is called by  RTreeSearch() to add selected node/line/area/isle to thelist */
+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 add selected node/line/area/isle to the list */
 int add_item(int id, struct ilist *list)
 {
     dig_list_add(list, id);
     return 1;
 }
 
+/* This function is called by RTreeSearch() to find an item in the list */
+int find_item(int id, struct ilist *list)
+{
+    dig_list_add(list, id);
+    return 0;
+}
+
 /*!
  * \brief Snap selected lines to existing vertex in threshold.
  *
@@ -117,6 +145,7 @@
     nvertices = 0;
     XPnts = NULL;
 
+    G_verbose_message(_("Snap vertices Pass 1: select points"));
     for (line_idx = 0; line_idx < List_lines->n_values; line_idx++) {
 	int v;
 
@@ -169,6 +198,9 @@
 
     /* 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;
@@ -236,6 +268,8 @@
 
     nsnapped = ncreated = 0;
 
+    G_verbose_message(_("Snap vertices Pass 3: snap to assigned points"));
+
     for (line_idx = 0; line_idx < List_lines->n_values; line_idx++) {
 	int v, spoint, anchor;
 	int changed = 0;
@@ -413,6 +447,9 @@
     G_free(Index);
     G_free(New);
     RTreeDestroyNode(RTree);
+
+    G_verbose_message(_("Snapped vertices: %d"), nsnapped);
+    G_verbose_message(_("New vertices: %d"), ncreated);
 }
 
 /* for qsort */
@@ -473,3 +510,397 @@
 
     return;
 }
+
+/*!
+   \brief Snap a line to reference lines in Map with threshold.
+  
+   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,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 *nsnapped, int *ncreated)
+{
+    struct line_pnts *LPoints, *NPoints;
+    struct line_cats *Cats;
+    int i, v, line, nlines;
+    int changed;
+    double thresh2;
+
+    int point;			/* index in points array */
+    int segment;		/* index in segments array */
+    int asegments;		/* number of allocated segments */
+    int apoints, nvertices;	/* number of allocated points and registered vertices */
+    XSEG *XSegs = NULL;		/* Array of segments */
+    XPNT *XPnts = NULL;		/* Array of points */
+    NEW2 *New = NULL;		/* Array of new points */
+    int anew = 0, nnew;		/* allocated new points , number of new points */
+    struct ilist *List;
+
+    struct Node *pnt_tree,	/* spatial index for reference points */
+                *seg_tree;	/* spatial index for reference segments */
+    struct Rect rect;
+
+    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();
+    pnt_tree = RTreeNewIndex();
+    seg_tree = RTreeNewIndex();
+
+    thresh2 = thresh * thresh;
+
+    point = segment = 1;	/* index starts from 1 ! */
+    nvertices = 0;
+    asegments = 0;
+    apoints = 0;
+
+    /* Add all vertices and all segments of all reference lines 
+     * to spatial indices */
+    nlines = reflist->n_values;
+    for (i = 0; i < nlines; i++) {
+
+	line = reflist->value[i];
+
+	G_debug(3, "line =  %d", line);
+	if (!Vect_line_alive(Map, line))
+	    continue;
+
+	Vect_read_line(Map, LPoints, Cats, line);
+	Vect_line_prune(LPoints);
+
+	for (v = 0; v < LPoints->n_points; v++) {
+	    G_debug(3, "  vertex v = %d", v);
+	    nvertices++;
+
+	    /* Box */
+	    rect.boundary[0] = LPoints->x[v];
+	    rect.boundary[3] = LPoints->x[v];
+	    rect.boundary[1] = LPoints->y[v];
+	    rect.boundary[4] = LPoints->y[v];
+
+	    /* Already registered ? */
+	    Vect_reset_list(List);
+	    RTreeSearch(pnt_tree, &rect, (void *)find_item, List);
+	    G_debug(3, "List : nvalues =  %d", List->n_values);
+
+	    if (List->n_values == 0) {	/* Not found */
+
+		/* Add to points tree */
+		RTreeInsertRect(&rect, point, &pnt_tree, 0);
+
+		if ((point - 1) == apoints) {
+		    apoints += 10000;
+		    XPnts =
+			(XPNT *) G_realloc(XPnts,
+					   (apoints + 1) * sizeof(XPNT));
+		}
+		XPnts[point].x = LPoints->x[v];
+		XPnts[point].y = LPoints->y[v];
+		XPnts[point].anchor = 0;
+
+		point++;
+	    }
+	    
+	    /* 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];
+		}
+
+		/* do not check for duplicates, too costly
+		 * because different segments can have identical boxes */
+		RTreeInsertRect(&rect, segment, &seg_tree, 0);
+
+		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];
+
+		segment++;
+	    }
+	}
+    }
+
+    /* go through all vertices of the line to snap */
+    for (v = 0; v < Points->n_points; v++) {
+	double dist2, tmpdist2;
+	double x, y;
+
+	dist2 = thresh2 + thresh2;
+	x = Points->x[v];
+	y = Points->y[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;
+
+	/* find nearest reference vertex */
+	Vect_reset_list(List);
+
+	RTreeSearch(pnt_tree, &rect, (void *)add_item, List);
+
+	for (i = 0; i < List->n_values; i++) {
+	    double dx, dy;
+	    
+	    point = List->value[i];
+	    
+	    dx = Points->x[v] - XPnts[point].x;
+	    dy = Points->y[v] - XPnts[point].y;
+	    
+	    tmpdist2 = dx * dx + dy * dy;
+	    
+	    if (tmpdist2 < dist2) {
+		dist2 = tmpdist2;
+
+		x = XPnts[point].x;
+		y = XPnts[point].y;
+	    }
+	}
+
+	if (dist2 <= thresh2 && dist2 > 0) {
+	    Points->x[v] = x;
+	    Points->y[v] = y;
+
+	    changed = 1;
+	    if (nsnapped)
+		(*nsnapped)++;
+	}
+    }
+
+    /* go through all vertices of the line to snap */
+    for (v = 0; v < Points->n_points; v++) {
+	double dist2, tmpdist2;
+	double x, y;
+	
+	dist2 = thresh2 + thresh2;
+	x = Points->x[v];
+	y = Points->y[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;
+
+	/* find nearest reference segment */
+	Vect_reset_list(List);
+
+	RTreeSearch(seg_tree, &rect, (void *)add_item, List);
+
+	for (i = 0; i < List->n_values; i++) {
+	    double tmpx, tmpy;
+	    int segment, status;
+	    
+	    segment = List->value[i];
+	    
+	    /* Check the distance */
+	    tmpdist2 =
+		dig_distance2_point_to_line(Points->x[v],
+					    Points->y[v],
+					    0, 
+					    XSegs[segment].x1,
+					    XSegs[segment].y1,
+					    0,
+					    XSegs[segment].x2,
+					    XSegs[segment].y2,
+					    0,
+					    0, &tmpx, &tmpy, NULL,
+					    NULL, &status);
+
+	    if (tmpdist2 < dist2 && status == 0) {
+		dist2 = tmpdist2;
+
+		x = tmpx;
+		y = tmpy;
+	    }
+	}
+
+	if (dist2 <= thresh2 && dist2 > 0) {
+	    Points->x[v] = x;
+	    Points->y[v] = y;
+	    
+	    changed = 1;
+	    if (nsnapped)
+		(*nsnapped)++;
+	}
+    }
+
+    RTreeDestroyNode(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;
+	double xmin, xmax, ymin, ymax;
+
+	x1 = Points->x[v];
+	x2 = Points->x[v + 1];
+	y1 = Points->y[v];
+	y2 = Points->y[v + 1];
+
+	Vect_append_point(NPoints, Points->x[v], Points->y[v],
+			  Points->z[v]);
+
+	/* 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;
+	}
+
+	rect.boundary[0] = xmin - thresh;
+	rect.boundary[3] = xmax + thresh;
+	rect.boundary[1] = ymin - thresh;
+	rect.boundary[4] = ymax + thresh;
+
+	/* Find points */
+	Vect_reset_list(List);
+	RTreeSearch(pnt_tree, &rect, (void *)add_item, List);
+
+	G_debug(3, "  %d points in box", List->n_values);
+
+	/* Snap to vertex in threshold different from end points */
+	nnew = 0;
+	for (i = 0; i < List->n_values; i++) {
+	    double dist2, along;
+	    int status;
+
+	    point = List->value[i];
+
+	    if (Points->x[v] == XPnts[i].x && 
+	        Points->y[v] == XPnts[i].y)
+		continue;	/* start point */
+
+	    if (Points->x[v + 1] == XPnts[i].x && 
+	        Points->y[v + 1] == XPnts[i].y)
+		continue;	/* end point */
+
+	    /* Check the distance */
+	    dist2 =
+		dig_distance2_point_to_line(XPnts[i].x,
+					    XPnts[i].y,
+					    0, x1, y1, 0,
+					    x2, y2, 0, 0, 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 = XPnts[i].x;
+		New[nnew].y = XPnts[i].y;
+		New[nnew].along = along;
+		nnew++;
+	    }
+	    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, 0);
+		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]);
+
+    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(LPoints);
+    Vect_destroy_line_struct(NPoints);
+    Vect_destroy_cats_struct(Cats);
+    Vect_destroy_list(List);
+    G_free(New);
+    G_free(XPnts);
+    RTreeDestroyNode(pnt_tree);
+
+    return changed;
+}



More information about the grass-commit mailing list