GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Vlib/snap.c
Go to the documentation of this file.
1 
20 #include <math.h>
21 #include <grass/gis.h>
22 #include <grass/Vect.h>
23 #include <grass/glocale.h>
24 
25 /* function prototypes */
26 static int sort_new(const void *pa, const void *pb);
27 
28 /* Vertex */
29 typedef struct
30 {
31  double x, y;
32  int anchor; /* 0 - anchor, do not snap this point, that means snap others to this */
33  /* >0 - index of anchor to which snap this point */
34  /* -1 - init value */
35 } XPNT;
36 
37 typedef struct
38 {
39  int anchor;
40  double along;
41 } NEW;
42 
43 /* This function is called by RTreeSearch() to add selected node/line/area/isle to thelist */
44 int add_item(int id, struct ilist *list)
45 {
46  dig_list_add(list, id);
47  return 1;
48 }
49 
68 /* As mentioned above, lines are not necessarily snapped to nearest vertex! For example:
69  |
70  | 1 line 3 is snapped to line 1,
71  | then line 2 is not snapped to common node at lines 1 and 3,
72  because it is already outside of threshold
73  ----------- 3
74 
75  |
76  | 2
77  |
78  */
79 void
80 Vect_snap_lines_list(struct Map_info *Map, struct ilist *List_lines,
81  double thresh, struct Map_info *Err)
82 {
83  struct line_pnts *Points, *NPoints;
84  struct line_cats *Cats;
85  int line, ltype, line_idx;
86  double thresh2;
87 
88  struct Node *RTree;
89  int point; /* index in points array */
90  int nanchors, ntosnap; /* number of anchors and number of points to be snapped */
91  int nsnapped, ncreated; /* number of snapped verices, number of new vertices (on segments) */
92  int apoints, npoints, nvertices; /* number of allocated points, registered points, vertices */
93  XPNT *XPnts; /* Array of points */
94  NEW *New = NULL; /* Array of new points */
95  int anew = 0, nnew; /* allocated new points , number of new points */
96  struct Rect rect;
97  struct ilist *List;
98  int *Index = NULL; /* indexes of anchors for vertices */
99  int aindex = 0; /* allocated Index */
100 
101  if (List_lines->n_values < 1)
102  return;
103 
104  Points = Vect_new_line_struct();
105  NPoints = Vect_new_line_struct();
106  Cats = Vect_new_cats_struct();
107  List = Vect_new_list();
108  RTree = RTreeNewIndex();
109 
110  thresh2 = thresh * thresh;
111 
112  /* Go through all lines in vector, and add each point to structure of points */
113  apoints = 0;
114  point = 1; /* index starts from 1 ! */
115  nvertices = 0;
116  XPnts = NULL;
117 
118  for (line_idx = 0; line_idx < List_lines->n_values; line_idx++) {
119  int v;
120 
121  G_percent(line_idx, List_lines->n_values, 2);
122 
123  line = List_lines->value[line_idx];
124 
125  G_debug(3, "line = %d", line);
126  if (!Vect_line_alive(Map, line))
127  continue;
128 
129  ltype = Vect_read_line(Map, Points, Cats, line);
130 
131  for (v = 0; v < Points->n_points; v++) {
132  G_debug(3, " vertex v = %d", v);
133  nvertices++;
134 
135  /* Box */
136  rect.boundary[0] = Points->x[v];
137  rect.boundary[3] = Points->x[v];
138  rect.boundary[1] = Points->y[v];
139  rect.boundary[4] = Points->y[v];
140  rect.boundary[2] = 0;
141  rect.boundary[5] = 0;
142 
143  /* Already registered ? */
144  Vect_reset_list(List);
145  RTreeSearch(RTree, &rect, (void *)add_item, List);
146  G_debug(3, "List : nvalues = %d", List->n_values);
147 
148  if (List->n_values == 0) { /* Not found */
149  /* Add to tree and to structure */
150  RTreeInsertRect(&rect, point, &RTree, 0);
151  if ((point - 1) == apoints) {
152  apoints += 10000;
153  XPnts =
154  (XPNT *) G_realloc(XPnts,
155  (apoints + 1) * sizeof(XPNT));
156  }
157  XPnts[point].x = Points->x[v];
158  XPnts[point].y = Points->y[v];
159  XPnts[point].anchor = -1;
160  point++;
161  }
162  }
163  }
164  G_percent(line_idx, List_lines->n_values, 2); /* finish it */
165 
166  npoints = point - 1;
167 
168  /* Go through all registered points and if not yet marked mark it as anchor and assign this anchor
169  * to all not yet marked points in threshold */
170  nanchors = ntosnap = 0;
171  for (point = 1; point <= npoints; point++) {
172  int i;
173 
174  G_percent(point, npoints, 2);
175 
176  G_debug(3, " point = %d", point);
177 
178  if (XPnts[point].anchor >= 0)
179  continue;
180 
181  XPnts[point].anchor = 0; /* make it anchor */
182  nanchors++;
183 
184  /* Find points in threshold */
185  rect.boundary[0] = XPnts[point].x - thresh;
186  rect.boundary[3] = XPnts[point].x + thresh;
187  rect.boundary[1] = XPnts[point].y - thresh;
188  rect.boundary[4] = XPnts[point].y + thresh;
189  rect.boundary[2] = 0;
190  rect.boundary[5] = 0;
191 
192  Vect_reset_list(List);
193  RTreeSearch(RTree, &rect, (void *)add_item, List);
194  G_debug(4, " %d points in threshold box", List->n_values);
195 
196  for (i = 0; i < List->n_values; i++) {
197  int pointb;
198  double dx, dy, dist2;
199 
200  pointb = List->value[i];
201  if (pointb == point)
202  continue;
203 
204  dx = XPnts[pointb].x - XPnts[point].x;
205  dy = XPnts[pointb].y - XPnts[point].y;
206  dist2 = dx * dx + dy * dy;
207 
208  if (dist2 > thresh2) /* outside threshold */
209  continue;
210 
211  /* doesn't have an anchor yet */
212  if (XPnts[pointb].anchor == -1) {
213  XPnts[pointb].anchor = point;
214  ntosnap++;
215  }
216  else if (XPnts[pointb].anchor > 0) { /* check distance to previously assigned anchor */
217  double dist2_a;
218 
219  dx = XPnts[XPnts[pointb].anchor].x - XPnts[pointb].x;
220  dy = XPnts[XPnts[pointb].anchor].y - XPnts[pointb].y;
221  dist2_a = dx * dx + dy * dy;
222 
223  /* replace old anchor */
224  if (dist2 < dist2_a) {
225  XPnts[pointb].anchor = point;
226  }
227  }
228  }
229  }
230 
231  /* Go through all lines and:
232  * 1) for all vertices: if not anchor snap it to its anchor
233  * 2) for all segments: snap it to all anchors in threshold (except anchors of vertices of course) */
234 
235  nsnapped = ncreated = 0;
236 
237  for (line_idx = 0; line_idx < List_lines->n_values; line_idx++) {
238  int v, spoint, anchor;
239  int changed = 0;
240 
241  G_percent(line_idx, List_lines->n_values, 2);
242 
243  line = List_lines->value[line_idx];
244 
245  G_debug(3, "line = %d", line);
246  if (!Vect_line_alive(Map, line))
247  continue;
248 
249  ltype = Vect_read_line(Map, Points, Cats, line);
250 
251  if (Points->n_points >= aindex) {
252  aindex = Points->n_points;
253  Index = (int *)G_realloc(Index, aindex * sizeof(int));
254  }
255 
256  /* Snap all vertices */
257  for (v = 0; v < Points->n_points; v++) {
258  /* Box */
259  rect.boundary[0] = Points->x[v];
260  rect.boundary[3] = Points->x[v];
261  rect.boundary[1] = Points->y[v];
262  rect.boundary[4] = Points->y[v];
263  rect.boundary[2] = 0;
264  rect.boundary[5] = 0;
265 
266  /* Find point ( should always find one point ) */
267  Vect_reset_list(List);
268 
269  RTreeSearch(RTree, &rect, (void *)add_item, List);
270 
271  spoint = List->value[0];
272  anchor = XPnts[spoint].anchor;
273 
274  if (anchor > 0) { /* to be snapped */
275  Points->x[v] = XPnts[anchor].x;
276  Points->y[v] = XPnts[anchor].y;
277  nsnapped++;
278  changed = 1;
279  Index[v] = anchor; /* point on new location */
280  }
281  else {
282  Index[v] = spoint; /* old point */
283  }
284  }
285 
286  /* New points */
287  Vect_reset_line(NPoints);
288 
289  /* Snap all segments to anchors in threshold */
290  for (v = 0; v < Points->n_points - 1; v++) {
291  int i;
292  double x1, x2, y1, y2, xmin, xmax, ymin, ymax;
293 
294  G_debug(3, " segment = %d end anchors : %d %d", v, Index[v],
295  Index[v + 1]);
296 
297  x1 = Points->x[v];
298  x2 = Points->x[v + 1];
299  y1 = Points->y[v];
300  y2 = Points->y[v + 1];
301 
302  Vect_append_point(NPoints, Points->x[v], Points->y[v],
303  Points->z[v]);
304 
305  /* Box */
306  if (x1 <= x2) {
307  xmin = x1;
308  xmax = x2;
309  }
310  else {
311  xmin = x2;
312  xmax = x1;
313  }
314  if (y1 <= y2) {
315  ymin = y1;
316  ymax = y2;
317  }
318  else {
319  ymin = y2;
320  ymax = y1;
321  }
322 
323  rect.boundary[0] = xmin - thresh;
324  rect.boundary[3] = xmax + thresh;
325  rect.boundary[1] = ymin - thresh;
326  rect.boundary[4] = ymax + thresh;
327  rect.boundary[2] = 0;
328  rect.boundary[5] = 0;
329 
330  /* Find points */
331  Vect_reset_list(List);
332  RTreeSearch(RTree, &rect, (void *)add_item, List);
333 
334  G_debug(3, " %d points in box", List->n_values);
335 
336  /* Snap to anchor in threshold different from end points */
337  nnew = 0;
338  for (i = 0; i < List->n_values; i++) {
339  double dist2, along;
340 
341  spoint = List->value[i];
342  G_debug(4, " spoint = %d anchor = %d", spoint,
343  XPnts[spoint].anchor);
344 
345  if (spoint == Index[v] || spoint == Index[v + 1])
346  continue; /* end point */
347  if (XPnts[spoint].anchor > 0)
348  continue; /* point is not anchor */
349 
350  /* Check the distance */
351  dist2 =
352  dig_distance2_point_to_line(XPnts[spoint].x,
353  XPnts[spoint].y, 0, x1, y1, 0,
354  x2, y2, 0, 0, NULL, NULL,
355  NULL, &along, NULL);
356 
357  G_debug(4, " distance = %lf", sqrt(dist2));
358 
359  if (dist2 <= thresh2) {
360  G_debug(4, " anchor in thresh, along = %lf", along);
361 
362  if (nnew == anew) {
363  anew += 100;
364  New = (NEW *) G_realloc(New, anew * sizeof(NEW));
365  }
366  New[nnew].anchor = spoint;
367  New[nnew].along = along;
368  nnew++;
369  }
370  }
371  G_debug(3, " nnew = %d", nnew);
372  /* insert new vertices */
373  if (nnew > 0) {
374  /* sort by distance along the segment */
375  qsort(New, sizeof(char) * nnew, sizeof(NEW), sort_new);
376 
377  for (i = 0; i < nnew; i++) {
378  anchor = New[i].anchor;
379  /* Vect_line_insert_point ( Points, ++v, XPnts[anchor].x, XPnts[anchor].y, 0); */
380  Vect_append_point(NPoints, XPnts[anchor].x,
381  XPnts[anchor].y, 0);
382  ncreated++;
383  }
384  changed = 1;
385  }
386  }
387 
388  /* append end point */
389  v = Points->n_points - 1;
390  Vect_append_point(NPoints, Points->x[v], Points->y[v], Points->z[v]);
391 
392  if (changed) { /* rewrite the line */
393  Vect_line_prune(NPoints); /* remove duplicates */
394  if (NPoints->n_points > 1 || ltype & GV_LINES) {
395  Vect_rewrite_line(Map, line, ltype, NPoints, Cats);
396  }
397  else {
398  Vect_delete_line(Map, line);
399  }
400  if (Err) {
401  Vect_write_line(Err, ltype, Points, Cats);
402  }
403  }
404  } /* for each line */
405  G_percent(line_idx, List_lines->n_values, 2); /* finish it */
406 
407  Vect_destroy_line_struct(Points);
408  Vect_destroy_line_struct(NPoints);
410  G_free(XPnts);
411  G_free(Index);
412  G_free(New);
413  RTreeDestroyNode(RTree);
414 }
415 
416 /* for qsort */
417 static int sort_new(const void *pa, const void *pb)
418 {
419  NEW *p1 = (NEW *) pa;
420  NEW *p2 = (NEW *) pb;
421 
422  if (p1->along < p2->along)
423  return -1;
424  if (p1->along > p2->along)
425  return 1;
426  return 1;
427 }
428 
441 void
442 Vect_snap_lines(struct Map_info *Map, int type, double thresh,
443  struct Map_info *Err)
444 {
445  int line, nlines, ltype;
446 
447  struct ilist *List;
448 
449  List = Vect_new_list();
450 
451  nlines = Vect_get_num_lines(Map);
452 
453  for (line = 1; line <= nlines; line++) {
454  G_debug(3, "line = %d", line);
455 
456  if (!Vect_line_alive(Map, line))
457  continue;
458 
459  ltype = Vect_read_line(Map, NULL, NULL, line);
460 
461  if (!(ltype & type))
462  continue;
463 
464  /* Vect_list_append(List, line); */
465  dig_list_add(List, line);
466  }
467 
468  Vect_snap_lines_list(Map, List, thresh, Err);
469 
470  Vect_destroy_list(List);
471 
472  return;
473 }