GRASS 8 Programmer's Manual 8.6.0dev(2026)-f6f2c534ea
Loading...
Searching...
No Matches
Vlib/snap.c
Go to the documentation of this file.
1/*!
2 * \file lib/vector/Vlib/snap.c
3 *
4 * \brief Vector library - Clean vector map (snap lines)
5 *
6 * Higher level functions for reading/writing/manipulating vectors.
7 *
8 * (C) 2001-2009 by the GRASS Development Team
9 *
10 * This program is free software under the GNU General Public License
11 * (>=v2). Read the file COPYING that comes with GRASS for details.
12 *
13 * \author Radim Blazek
14 * \author update to GRASS 7 Markus Metz
15 */
16
17#include <errno.h>
18#include <stdlib.h>
19#include <string.h>
20#include <sys/stat.h>
21#include <fcntl.h>
22#include <unistd.h>
23#include <math.h>
24#include <grass/vector.h>
25#include <grass/gis.h>
26#include <grass/glocale.h>
27#include <grass/kdtree.h>
28
29/* translate segment to box and back */
30#define X1W 0x01 /* x1 is West, x2 East */
31#define Y1S 0x02 /* y1 is South, y2 North */
32#define Z1B 0x04 /* z1 is Bottom, z2 Top */
33
34/* Vertex */
35typedef struct {
36 double x, y, z;
37 int anchor; /* 0 - anchor, do not snap this point, that means snap others to
38 this */
39 /* >0 - index of anchor to which snap this point */
40 /* -1 - init value */
41} XPNT;
42
43typedef struct {
44 int anchor;
45 double along;
46} NEW;
47
48/* for qsort */
49static int sort_new(const void *pa, const void *pb)
50{
51 NEW *p1 = (NEW *)pa;
52 NEW *p2 = (NEW *)pb;
53
54 return (p1->along < p2->along ? -1 : (p1->along > p2->along));
55
56 /*
57 if (p1->along < p2->along)
58 return -1;
59 if (p1->along > p2->along)
60 return 1;
61 return 1;
62 */
63}
64
65typedef struct {
66 double x, y, z, along;
67} NEW2;
68
69/* for qsort */
70static int sort_new2(const void *pa, const void *pb)
71{
72 NEW2 *p1 = (NEW2 *)pa;
73 NEW2 *p2 = (NEW2 *)pb;
74
75 return (p1->along < p2->along ? -1 : (p1->along > p2->along));
76}
77
78/* This function is called by RTreeSearch() to find a vertex */
79static int find_item(int id, const struct RTree_Rect *rect UNUSED, void *list)
80{
81 G_ilist_add((struct ilist *)list, id);
82 return 0;
83}
84
85/* This function is called by RTreeSearch() to add selected node/line/area/isle
86 * to the list */
87static int add_item(int id, const struct RTree_Rect *rect UNUSED, void *list)
88{
89 G_ilist_add((struct ilist *)list, id);
90 return 1;
91}
92
93/* This function is called by RTreeSearch() to add selected node/line/area/isle
94 * to the list */
95static int find_item_box(int id, const struct RTree_Rect *rect, void *list)
96{
97 struct bound_box box;
98
99 box.W = rect->boundary[0];
100 box.S = rect->boundary[1];
101 box.B = rect->boundary[2];
102 box.E = rect->boundary[3];
103 box.N = rect->boundary[4];
104 box.T = rect->boundary[5];
105
106 dig_boxlist_add((struct boxlist *)list, id, &box);
107
108 return 0;
109}
110
111/* This function is called by RTreeSearch() to add selected node/line/area/isle
112 * to the list */
113static int add_item_box(int id, const struct RTree_Rect *rect, void *list)
114{
115 struct bound_box box;
116
117 box.W = rect->boundary[0];
118 box.S = rect->boundary[1];
119 box.B = rect->boundary[2];
120 box.E = rect->boundary[3];
121 box.N = rect->boundary[4];
122 box.T = rect->boundary[5];
123
124 dig_boxlist_add((struct boxlist *)list, id, &box);
125
126 return 1;
127}
128
129static void Vect_snap_lines_list_rtree(struct Map_info *, const struct ilist *,
130 double, struct Map_info *);
131
132static void Vect_snap_lines_list_kdtree(struct Map_info *, const struct ilist *,
133 double, struct Map_info *);
134
135/*!
136 \brief Snap selected lines to existing vertex in threshold.
137
138 Snap selected lines to existing vertices of other selected lines.
139 3D snapping is not supported.
140
141 Lines showing how vertices were snapped may be optionally written to error
142 map. Input map must be opened on level 2 for update at least on
143 GV_BUILD_BASE.
144
145 As mentioned above, lines are not necessarily snapped to nearest vertex! For
146 example: <pre>
147 |
148 | 1 line 3 is snapped to line 1,
149 | then line 2 is not snapped to common node at lines 1 and 3,
150 because it is already outside of threshold
151 ----------- 3
152
153 |
154 | 2
155 |
156 </pre>
157
158 The algorithm selects anchor vertices and snaps non-anchor vertices
159 to these anchors.
160 The distance between anchor vertices is always > threshold.
161 If there is more than one anchor vertex within threshold around a
162 non-anchor vertex, this vertex is snapped to the nearest anchor
163 vertex within threshold.
164
165 \param Map input map where vertices will be snapped
166 \param List_lines list of lines to snap
167 \param thresh threshold in which snap vertices
168 \param[out] Err vector map where lines representing snap are written or NULL
169
170 \return void
171 */
173 double thresh, struct Map_info *Err)
174{
175 if (getenv("GRASS_VECTOR_LOWMEM"))
176 Vect_snap_lines_list_rtree(Map, List_lines, thresh, Err);
177 else
178 Vect_snap_lines_list_kdtree(Map, List_lines, thresh, Err);
179}
180
181static void Vect_snap_lines_list_kdtree(struct Map_info *Map,
182 const struct ilist *List_lines,
183 double thresh, struct Map_info *Err)
184{
185 struct line_pnts *Points, *NPoints;
186 struct line_cats *Cats;
187 int line, ltype, line_idx;
188 double thresh2;
189
190 int point; /* index in points array */
191 int nsnapped, ncreated; /* number of snapped verices, number of new vertices
192 (on segments) */
193 int apoints, npoints; /* number of allocated points, registered points */
194 XPNT *XPnts; /* Array of points */
195 NEW *New = NULL; /* Array of new points */
196 int anew = 0, nnew; /* allocated new points, number of new points */
197 struct ilist *List;
198 int *Index = NULL; /* indexes of anchors for vertices */
199 int aindex = 0; /* allocated Index */
200
201 struct kdtree *KDTree;
202 double c[2];
203 double *kdd;
204 int *kduid, kd_found;
205
206 if (List_lines->n_values < 1)
207 return;
208
209 Points = Vect_new_line_struct();
213
215
217
218 /* Go through all lines in vector, and add each point to structure of points
219 */
220 apoints = 0;
221 point = 1; /* index starts from 1 ! */
222 XPnts = NULL;
223
224 G_important_message(_("Snap vertices Pass 1: select points"));
225 for (line_idx = 0; line_idx < List_lines->n_values; line_idx++) {
226 int v;
227
228 G_percent(line_idx, List_lines->n_values, 2);
229
230 line = List_lines->value[line_idx];
231
232 G_debug(3, "line = %d", line);
233 if (!Vect_line_alive(Map, line))
234 continue;
235
236 ltype = Vect_read_line(Map, Points, Cats, line);
237
238 for (v = 0; v < Points->n_points; v++) {
239
240 G_debug(3, " vertex v = %d", v);
241
242 /* coords */
243 c[0] = Points->x[v];
244 c[1] = Points->y[v];
245
246 if (kdtree_insert(KDTree, c, point, 0)) {
247 /* Add to structure */
248 if ((point - 1) == apoints) {
249 apoints += 10000;
250 XPnts =
251 (XPNT *)G_realloc(XPnts, (apoints + 1) * sizeof(XPNT));
252 }
253 XPnts[point].x = Points->x[v];
254 XPnts[point].y = Points->y[v];
255 XPnts[point].anchor = -1;
256 point++;
257 }
258 }
259 }
260 G_percent(line_idx, List_lines->n_values, 2); /* finish it */
261
262 npoints = point - 1;
263
264 /* Go through all registered points and if not yet marked mark it as anchor
265 * and assign this anchor to all not yet marked points in threshold */
266
267 G_important_message(_("Snap vertices Pass 2: assign anchor vertices"));
268
269 for (point = 1; point <= npoints; point++) {
270 int i;
271
272 G_percent(point, npoints, 4);
273
274 G_debug(3, " point = %d", point);
275
276 if (XPnts[point].anchor >= 0)
277 continue;
278
279 XPnts[point].anchor = 0; /* make it anchor */
280
281 /* Find points in threshold */
282 c[0] = XPnts[point].x;
283 c[1] = XPnts[point].y;
284
286 kd_found = kdtree_dnn(KDTree, c, &kduid, &kdd, thresh, &point);
287 G_debug(4, " %d points in threshold box", kd_found);
288
289 for (i = 0; i < kd_found; i++) {
290 int pointb;
291 double dx, dy, dist2;
292
293 pointb = kduid[i];
294 if (pointb == point)
295 continue;
296
297 dx = XPnts[pointb].x - XPnts[point].x;
298 dy = XPnts[pointb].y - XPnts[point].y;
299 dist2 = dx * dx + dy * dy;
300
301 if (dist2 > thresh2) /* outside threshold */
302 continue;
303
304 /* doesn't have an anchor yet */
305 if (XPnts[pointb].anchor == -1) {
306 XPnts[pointb].anchor = point;
307 }
308 else if (XPnts[pointb].anchor >
309 0) { /* check distance to previously assigned anchor */
310 double dist2_a;
311
312 dx = XPnts[XPnts[pointb].anchor].x - XPnts[pointb].x;
313 dy = XPnts[XPnts[pointb].anchor].y - XPnts[pointb].y;
314 dist2_a = dx * dx + dy * dy;
315
316 /* replace old anchor */
317 if (dist2 < dist2_a) {
318 XPnts[pointb].anchor = point;
319 }
320 }
321 }
322 if (kd_found) {
323 G_free(kdd);
324 G_free(kduid);
325 }
326 }
327
328 /* Go through all lines and:
329 * 1) for all vertices: if not anchor snap it to its anchor
330 * 2) for all segments: snap it to all anchors in threshold (except
331 * anchors of vertices of course) */
332
333 nsnapped = ncreated = 0;
334
335 G_important_message(_("Snap vertices Pass 3: snap to assigned points"));
336
337 for (line_idx = 0; line_idx < List_lines->n_values; line_idx++) {
338 int v, spoint, anchor;
339 int changed = 0;
340 double kddist;
341
342 G_percent(line_idx, List_lines->n_values, 2);
343
344 line = List_lines->value[line_idx];
345
346 G_debug(3, "line = %d", line);
347 if (!Vect_line_alive(Map, line))
348 continue;
349
350 ltype = Vect_read_line(Map, Points, Cats, line);
351
352 if (Points->n_points >= aindex) {
353 aindex = Points->n_points;
354 Index = (int *)G_realloc(Index, aindex * sizeof(int));
355 }
356
357 /* Snap all vertices */
358 G_debug(3, "Snap all vertices");
359 for (v = 0; v < Points->n_points; v++) {
360 /* Box */
361 c[0] = Points->x[v];
362 c[1] = Points->y[v];
363
364 /* Find point ( should always find one point ) */
366
367 spoint = -1;
368 kdtree_knn(KDTree, c, &spoint, &kddist, 1, NULL);
369 if (spoint == -1)
370 G_fatal_error("Point not in KD Tree");
371
372 anchor = XPnts[spoint].anchor;
373
374 if (anchor > 0) { /* to be snapped */
375 Points->x[v] = XPnts[anchor].x;
376 Points->y[v] = XPnts[anchor].y;
377 nsnapped++;
378 changed = 1;
379 Index[v] = anchor; /* point on new location */
380 }
381 else {
382 Index[v] = spoint; /* old point */
383 }
384 }
385
386 /* New points */
388
389 /* Snap all segments to anchors in threshold */
390 G_debug(3, "Snap all segments");
391 for (v = 0; v < Points->n_points - 1; v++) {
392 int i;
393 double x1, x2, y1, y2, xmin, xmax, ymin, ymax;
394 double rc[4];
395
396 G_debug(3, " segment = %d end anchors : %d %d", v, Index[v],
397 Index[v + 1]);
398
399 x1 = Points->x[v];
400 x2 = Points->x[v + 1];
401 y1 = Points->y[v];
402 y2 = Points->y[v + 1];
403
404 Vect_append_point(NPoints, Points->x[v], Points->y[v],
405 Points->z[v]);
406
407 /* Box */
408 if (x1 <= x2) {
409 xmin = x1;
410 xmax = x2;
411 }
412 else {
413 xmin = x2;
414 xmax = x1;
415 }
416 if (y1 <= y2) {
417 ymin = y1;
418 ymax = y2;
419 }
420 else {
421 ymin = y2;
422 ymax = y1;
423 }
424
425 /* Find points */
427 G_debug(3, " search anchors for segment %g,%g to %g,%g", x1, y1,
428 x2, y2);
429 /* distance search: circle around midpoint encompassing
430 * endpoints
431 * box search: box encompassing endpoints,
432 * smaller than corresponding circle */
433 rc[0] = xmin - thresh * 2;
434 rc[1] = ymin - thresh * 2;
435 rc[2] = xmax + thresh * 2;
436 rc[3] = ymax + thresh * 2;
437
439
440 G_debug(3, " %d points in box", kd_found);
441
442 /* Snap to anchor in threshold different from end points */
443 nnew = 0;
444 for (i = 0; i < kd_found; i++) {
445 double dist2, along;
446 int status;
447
448 spoint = kduid[i];
449 G_debug(4, " spoint = %d anchor = %d", spoint,
450 XPnts[spoint].anchor);
451
452 if (spoint == Index[v] || spoint == Index[v + 1])
453 continue; /* end point */
454 if (XPnts[spoint].anchor > 0)
455 continue; /* point is not anchor */
456
457 /* Check the distance */
459 XPnts[spoint].x, XPnts[spoint].y, 0, x1, y1, 0, x2, y2, 0,
460 0, NULL, NULL, NULL, &along, &status);
461
462 G_debug(4, " distance = %lf", sqrt(dist2));
463
464 if (status == 0 && dist2 <= thresh2) {
465 G_debug(4, " anchor in thresh, along = %lf", along);
466
467 if (nnew == anew) {
468 anew += 100;
469 New = (NEW *)G_realloc(New, anew * sizeof(NEW));
470 }
471 New[nnew].anchor = spoint;
472 New[nnew].along = along;
473 nnew++;
474 }
475 }
476 if (kd_found) {
477 G_free(kduid);
478 }
479 G_debug(3, " nnew = %d", nnew);
480 /* insert new vertices */
481 if (nnew > 0) {
482 /* sort by distance along the segment */
483 qsort(New, sizeof(char) * nnew, sizeof(NEW), sort_new);
484
485 for (i = 0; i < nnew; i++) {
486 anchor = New[i].anchor;
487 /* Vect_line_insert_point ( Points, ++v, XPnts[anchor].x,
488 * XPnts[anchor].y, 0); */
489 Vect_append_point(NPoints, XPnts[anchor].x, XPnts[anchor].y,
490 0);
491 ncreated++;
492 }
493 changed = 1;
494 }
495 }
496
497 /* append end point */
498 v = Points->n_points - 1;
499 Vect_append_point(NPoints, Points->x[v], Points->y[v], Points->z[v]);
500
501 if (changed) { /* rewrite the line */
502 Vect_line_prune(NPoints); /* remove duplicates */
503 if (NPoints->n_points > 1 || !(ltype & GV_LINES)) {
505 }
506 else {
507 Vect_delete_line(Map, line);
508 }
509 if (Err) {
510 Vect_write_line(Err, ltype, Points, Cats);
511 }
512 }
513 } /* for each line */
514 G_percent(line_idx, List_lines->n_values, 2); /* finish it */
515
519 G_free(XPnts);
520 G_free(Index);
521 G_free(New);
523
524 G_verbose_message(_("Snapped vertices: %d"), nsnapped);
525 G_verbose_message(_("New vertices: %d"), ncreated);
527}
528
529static void Vect_snap_lines_list_rtree(struct Map_info *Map,
530 const struct ilist *List_lines,
531 double thresh, struct Map_info *Err)
532{
533 struct line_pnts *Points, *NPoints;
534 struct line_cats *Cats;
535 int line, ltype, line_idx;
536 double thresh2;
537
538 int point; /* index in points array */
539 int nsnapped, ncreated; /* number of snapped verices, number of new vertices
540 (on segments) */
541 int apoints, npoints; /* number of allocated points, registered points */
542 XPNT *XPnts; /* Array of points */
543 NEW *New = NULL; /* Array of new points */
544 int anew = 0, nnew; /* allocated new points , number of new points */
545 struct ilist *List;
546 int *Index = NULL; /* indexes of anchors for vertices */
547 int aindex = 0; /* allocated Index */
548
549 struct RTree *RTree;
550 int rtreefd = -1;
551 static struct RTree_Rect rect;
552 static int rect_init = 0;
553
554 if (!rect_init) {
555 rect.boundary = G_malloc(6 * sizeof(RectReal));
556 rect_init = 6;
557 }
558
559 if (List_lines->n_values < 1)
560 return;
561
562 Points = Vect_new_line_struct();
566 if (getenv("GRASS_VECTOR_LOWMEM")) {
567 char *filename = G_tempfile();
568
569 rtreefd = open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
570 if (rtreefd < 0) {
571 G_fatal_error(_("Unable to create temporary file <%s>: %s"),
572 filename, strerror(errno));
573 }
574 if (remove(filename) != 0) {
575 G_warning(_("Unable to remove temporary file <%s>: %s"), filename,
576 strerror(errno));
577 }
578 G_free(filename);
579 }
581
583
584 /* Go through all lines in vector, and add each point to structure of points
585 */
586 apoints = 0;
587 point = 1; /* index starts from 1 ! */
588 XPnts = NULL;
589
590 G_important_message(_("Snap vertices Pass 1: select points"));
591 for (line_idx = 0; line_idx < List_lines->n_values; line_idx++) {
592 int v;
593
594 G_percent(line_idx, List_lines->n_values, 2);
595
596 line = List_lines->value[line_idx];
597
598 G_debug(3, "line = %d", line);
599 if (!Vect_line_alive(Map, line))
600 continue;
601
602 ltype = Vect_read_line(Map, Points, Cats, line);
603
604 for (v = 0; v < Points->n_points; v++) {
605 G_debug(3, " vertex v = %d", v);
606
607 /* Box */
608 rect.boundary[0] = Points->x[v];
609 rect.boundary[3] = Points->x[v];
610 rect.boundary[1] = Points->y[v];
611 rect.boundary[4] = Points->y[v];
612 rect.boundary[2] = 0;
613 rect.boundary[5] = 0;
614
615 /* Already registered ? */
617 RTreeSearch(RTree, &rect, find_item, List);
618 G_debug(3, "List : nvalues = %d", List->n_values);
619
620 if (List->n_values == 0) { /* Not found */
621 /* Add to tree and to structure */
622 RTreeInsertRect(&rect, point, RTree);
623 if ((point - 1) == apoints) {
624 apoints += 10000;
625 XPnts =
626 (XPNT *)G_realloc(XPnts, (apoints + 1) * sizeof(XPNT));
627 }
628 XPnts[point].x = Points->x[v];
629 XPnts[point].y = Points->y[v];
630 XPnts[point].anchor = -1;
631 point++;
632 }
633 }
634 }
635 G_percent(line_idx, List_lines->n_values, 2); /* finish it */
636
637 npoints = point - 1;
638
639 /* Go through all registered points and if not yet marked mark it as anchor
640 * and assign this anchor to all not yet marked points in threshold */
641
642 G_important_message(_("Snap vertices Pass 2: assign anchor vertices"));
643
644 for (point = 1; point <= npoints; point++) {
645 int i;
646
647 G_percent(point, npoints, 4);
648
649 G_debug(3, " point = %d", point);
650
651 if (XPnts[point].anchor >= 0)
652 continue;
653
654 XPnts[point].anchor = 0; /* make it anchor */
655
656 /* Find points in threshold */
657 rect.boundary[0] = XPnts[point].x - thresh;
658 rect.boundary[3] = XPnts[point].x + thresh;
659 rect.boundary[1] = XPnts[point].y - thresh;
660 rect.boundary[4] = XPnts[point].y + thresh;
661 rect.boundary[2] = 0;
662 rect.boundary[5] = 0;
663
665 RTreeSearch(RTree, &rect, add_item, List);
666 G_debug(4, " %d points in threshold box", List->n_values);
667
668 for (i = 0; i < List->n_values; i++) {
669 int pointb;
670 double dx, dy, dist2;
671
672 pointb = List->value[i];
673 if (pointb == point)
674 continue;
675
676 dx = XPnts[pointb].x - XPnts[point].x;
677 dy = XPnts[pointb].y - XPnts[point].y;
678 dist2 = dx * dx + dy * dy;
679
680 if (dist2 > thresh2) /* outside threshold */
681 continue;
682
683 /* doesn't have an anchor yet */
684 if (XPnts[pointb].anchor == -1) {
685 XPnts[pointb].anchor = point;
686 }
687 else if (XPnts[pointb].anchor >
688 0) { /* check distance to previously assigned anchor */
689 double dist2_a;
690
691 dx = XPnts[XPnts[pointb].anchor].x - XPnts[pointb].x;
692 dy = XPnts[XPnts[pointb].anchor].y - XPnts[pointb].y;
693 dist2_a = dx * dx + dy * dy;
694
695 /* replace old anchor */
696 if (dist2 < dist2_a) {
697 XPnts[pointb].anchor = point;
698 }
699 }
700 }
701 }
702
703 /* Go through all lines and:
704 * 1) for all vertices: if not anchor snap it to its anchor
705 * 2) for all segments: snap it to all anchors in threshold (except
706 * anchors of vertices of course) */
707
708 nsnapped = ncreated = 0;
709
710 G_important_message(_("Snap vertices Pass 3: snap to assigned points"));
711
712 for (line_idx = 0; line_idx < List_lines->n_values; line_idx++) {
713 int v, spoint, anchor;
714 int changed = 0;
715
716 G_percent(line_idx, List_lines->n_values, 2);
717
718 line = List_lines->value[line_idx];
719
720 G_debug(3, "line = %d", line);
721 if (!Vect_line_alive(Map, line))
722 continue;
723
724 ltype = Vect_read_line(Map, Points, Cats, line);
725
726 if (Points->n_points >= aindex) {
727 aindex = Points->n_points;
728 Index = (int *)G_realloc(Index, aindex * sizeof(int));
729 }
730
731 /* Snap all vertices */
732 for (v = 0; v < Points->n_points; v++) {
733 /* Box */
734 rect.boundary[0] = Points->x[v];
735 rect.boundary[3] = Points->x[v];
736 rect.boundary[1] = Points->y[v];
737 rect.boundary[4] = Points->y[v];
738 rect.boundary[2] = 0;
739 rect.boundary[5] = 0;
740
741 /* Find point ( should always find one point ) */
743
744 RTreeSearch(RTree, &rect, add_item, List);
745
746 spoint = List->value[0];
747 anchor = XPnts[spoint].anchor;
748
749 if (anchor > 0) { /* to be snapped */
750 Points->x[v] = XPnts[anchor].x;
751 Points->y[v] = XPnts[anchor].y;
752 nsnapped++;
753 changed = 1;
754 Index[v] = anchor; /* point on new location */
755 }
756 else {
757 Index[v] = spoint; /* old point */
758 }
759 }
760
761 /* New points */
763
764 /* Snap all segments to anchors in threshold */
765 for (v = 0; v < Points->n_points - 1; v++) {
766 int i;
767 double x1, x2, y1, y2, xmin, xmax, ymin, ymax;
768
769 G_debug(3, " segment = %d end anchors : %d %d", v, Index[v],
770 Index[v + 1]);
771
772 x1 = Points->x[v];
773 x2 = Points->x[v + 1];
774 y1 = Points->y[v];
775 y2 = Points->y[v + 1];
776
777 Vect_append_point(NPoints, Points->x[v], Points->y[v],
778 Points->z[v]);
779
780 /* Box */
781 if (x1 <= x2) {
782 xmin = x1;
783 xmax = x2;
784 }
785 else {
786 xmin = x2;
787 xmax = x1;
788 }
789 if (y1 <= y2) {
790 ymin = y1;
791 ymax = y2;
792 }
793 else {
794 ymin = y2;
795 ymax = y1;
796 }
797
798 rect.boundary[0] = xmin - thresh;
799 rect.boundary[3] = xmax + thresh;
800 rect.boundary[1] = ymin - thresh;
801 rect.boundary[4] = ymax + thresh;
802 rect.boundary[2] = 0;
803 rect.boundary[5] = 0;
804
805 /* Find points */
807 RTreeSearch(RTree, &rect, add_item, List);
808
809 G_debug(3, " %d points in box", List->n_values);
810
811 /* Snap to anchor in threshold different from end points */
812 nnew = 0;
813 for (i = 0; i < List->n_values; i++) {
814 double dist2, along;
815 int status;
816
817 spoint = List->value[i];
818 G_debug(4, " spoint = %d anchor = %d", spoint,
819 XPnts[spoint].anchor);
820
821 if (spoint == Index[v] || spoint == Index[v + 1])
822 continue; /* end point */
823 if (XPnts[spoint].anchor > 0)
824 continue; /* point is not anchor */
825
826 /* Check the distance */
828 XPnts[spoint].x, XPnts[spoint].y, 0, x1, y1, 0, x2, y2, 0,
829 0, NULL, NULL, NULL, &along, &status);
830
831 G_debug(4, " distance = %lf", sqrt(dist2));
832
833 if (status == 0 && dist2 <= thresh2) {
834 G_debug(4, " anchor in thresh, along = %lf", along);
835
836 if (nnew == anew) {
837 anew += 100;
838 New = (NEW *)G_realloc(New, anew * sizeof(NEW));
839 }
840 New[nnew].anchor = spoint;
841 New[nnew].along = along;
842 nnew++;
843 }
844 }
845 G_debug(3, " nnew = %d", nnew);
846 /* insert new vertices */
847 if (nnew > 0) {
848 /* sort by distance along the segment */
849 qsort(New, sizeof(char) * nnew, sizeof(NEW), sort_new);
850
851 for (i = 0; i < nnew; i++) {
852 anchor = New[i].anchor;
853 /* Vect_line_insert_point ( Points, ++v, XPnts[anchor].x,
854 * XPnts[anchor].y, 0); */
855 Vect_append_point(NPoints, XPnts[anchor].x, XPnts[anchor].y,
856 0);
857 ncreated++;
858 }
859 changed = 1;
860 }
861 }
862
863 /* append end point */
864 v = Points->n_points - 1;
865 Vect_append_point(NPoints, Points->x[v], Points->y[v], Points->z[v]);
866
867 if (changed) { /* rewrite the line */
868 Vect_line_prune(NPoints); /* remove duplicates */
869 if (NPoints->n_points > 1 || !(ltype & GV_LINES)) {
871 }
872 else {
873 Vect_delete_line(Map, line);
874 }
875 if (Err) {
876 Vect_write_line(Err, ltype, Points, Cats);
877 }
878 }
879 } /* for each line */
880 G_percent(line_idx, List_lines->n_values, 2); /* finish it */
881
885 G_free(XPnts);
886 G_free(Index);
887 G_free(New);
889 if (rtreefd >= 0)
890 close(rtreefd);
891
892 G_verbose_message(_("Snapped vertices: %d"), nsnapped);
893 G_verbose_message(_("New vertices: %d"), ncreated);
895}
896
897/*!
898 \brief Snap lines in vector map to existing vertex in threshold.
899
900 For details see Vect_snap_lines_list()
901
902 \param[in] Map input map where vertices will be snapped
903 \param[in] type type of lines to snap
904 \param[in] thresh threshold in which snap vertices
905 \param[out] Err vector map where lines representing snap are written or NULL
906
907 \return void
908 */
909void Vect_snap_lines(struct Map_info *Map, int type, double thresh,
910 struct Map_info *Err)
911{
912 int line, nlines, ltype;
913 struct ilist *List;
914
916
917 nlines = Vect_get_num_lines(Map);
918
919 G_important_message(_("Reading features..."));
920 for (line = 1; line <= nlines; line++) {
921 G_debug(3, "line = %d", line);
922
923 if (!Vect_line_alive(Map, line))
924 continue;
925
926 ltype = Vect_read_line(Map, NULL, NULL, line);
927
928 if (!(ltype & type))
929 continue;
930
931 G_ilist_add(List, line);
932 }
933
935
937
938 return;
939}
940
941/*!
942 \brief Snap a line to reference lines in Map with threshold.
943
944 3D snapping is supported. The line to snap and the reference lines
945 can but do not need to be in different vector maps.
946
947 Vect_snap_line() uses less memory, but is slower than
948 Vect_snap_lines_list()
949
950 For details on snapping, see Vect_snap_lines_list()
951
952 \param[in] Map input map with reference lines
953 \param[in] reflist list of reference lines
954 \param[in,out] Points line points to snap
955 \param[in] thresh threshold in which to snap vertices
956 \param[in] with_z 2D or 3D snapping
957 \param[in,out] nsnapped number of snapped vertices
958 \param[in,out] ncreated number of new vertices (on segments)
959
960 \return 1 if line was changed, otherwise 0
961 */
963 struct line_pnts *Points, double thresh, int with_z,
964 int *nsnapped, int *ncreated)
965{
966 struct line_pnts *LPoints, *NPoints;
967 struct line_cats *Cats;
968 int i, v, line, nlines;
969 int changed;
970 double thresh2;
971
972 int point; /* index in points array */
973 int segment; /* index in segments array */
974 int asegments; /* number of allocated segments */
975 char *XSegs = NULL; /* Array of segments */
976 NEW2 *New = NULL; /* Array of new points */
977 int anew = 0, nnew; /* allocated new points , number of new points */
978 struct boxlist *List;
979
980 struct RTree *pnt_tree, /* spatial index for reference points */
981 *seg_tree; /* spatial index for reference segments */
982 struct RTree_Rect rect;
983
984 rect.boundary = G_malloc(6 * sizeof(RectReal));
985 rect.boundary[0] = 0;
986 rect.boundary[1] = 0;
987 rect.boundary[2] = 0;
988 rect.boundary[3] = 0;
989 rect.boundary[4] = 0;
990 rect.boundary[5] = 0;
991
992 changed = 0;
993 if (nsnapped)
994 *nsnapped = 0;
995 if (ncreated)
996 *ncreated = 0;
997
998 point = Points->n_points;
999 Vect_line_prune(Points);
1000 if (point != Points->n_points)
1001 changed = 1;
1002
1003 nlines = reflist->n_values;
1004 if (nlines < 1) {
1005 G_free(rect.boundary);
1006 return changed;
1007 }
1008
1013 with_z = (with_z != 0);
1014 pnt_tree = RTreeCreateTree(-1, 0, 2 + with_z);
1016 seg_tree = RTreeCreateTree(-1, 0, 2 + with_z);
1018
1019 thresh2 = thresh * thresh;
1020
1021 point = segment = 1; /* index starts from 1 ! */
1022 asegments = 0;
1023
1024 /* Add all vertices and all segments of all reference lines
1025 * to spatial indices */
1026 nlines = reflist->n_values;
1027 for (i = 0; i < nlines; i++) {
1028
1029 line = reflist->value[i];
1030
1031 G_debug(3, "line = %d", line);
1032 if (!Vect_line_alive(Map, line))
1033 continue;
1034
1035 Vect_read_line(Map, LPoints, Cats, line);
1037
1038 for (v = 0; v < LPoints->n_points; v++) {
1039 G_debug(3, " vertex v = %d", v);
1040
1041 /* Box */
1042 rect.boundary[0] = LPoints->x[v];
1043 rect.boundary[3] = LPoints->x[v];
1044 rect.boundary[1] = LPoints->y[v];
1045 rect.boundary[4] = LPoints->y[v];
1046 if (with_z) {
1047 rect.boundary[2] = LPoints->z[v];
1048 rect.boundary[5] = LPoints->z[v];
1049 }
1050
1051 /* Already registered ? */
1053 RTreeSearch(pnt_tree, &rect, find_item_box, (void *)List);
1054 G_debug(3, "List : nvalues = %d", List->n_values);
1055
1056 if (List->n_values == 0) { /* Not found */
1057
1058 /* Add to points tree */
1059 RTreeInsertRect(&rect, point, pnt_tree);
1060
1061 point++;
1062 }
1063
1064 /* reference segments */
1065 if (v) {
1066 char sides = 0;
1067
1068 /* Box */
1069 if (LPoints->x[v - 1] < LPoints->x[v]) {
1070 rect.boundary[0] = LPoints->x[v - 1];
1071 rect.boundary[3] = LPoints->x[v];
1072 sides |= X1W;
1073 }
1074 else {
1075 rect.boundary[0] = LPoints->x[v];
1076 rect.boundary[3] = LPoints->x[v - 1];
1077 }
1078 if (LPoints->y[v - 1] < LPoints->y[v]) {
1079 rect.boundary[1] = LPoints->y[v - 1];
1080 rect.boundary[4] = LPoints->y[v];
1081 sides |= Y1S;
1082 }
1083 else {
1084 rect.boundary[1] = LPoints->y[v];
1085 rect.boundary[4] = LPoints->y[v - 1];
1086 }
1087 if (LPoints->z[v - 1] < LPoints->z[v]) {
1088 rect.boundary[2] = LPoints->z[v - 1];
1089 rect.boundary[5] = LPoints->z[v];
1090 sides |= Z1B;
1091 }
1092 else {
1093 rect.boundary[2] = LPoints->z[v];
1094 rect.boundary[5] = LPoints->z[v - 1];
1095 }
1096
1097 /* do not check for duplicates, too costly
1098 * because different segments can have identical boxes */
1099 RTreeInsertRect(&rect, segment, seg_tree);
1100
1101 if ((segment - 1) == asegments) {
1102 asegments += 1000;
1103 XSegs = (char *)G_realloc(XSegs,
1104 (asegments + 1) * sizeof(char));
1105 }
1106 XSegs[segment] = sides;
1107 segment++;
1108 }
1109 }
1110 }
1111
1112 /* go through all vertices of the line to snap */
1113 /* find nearest reference vertex */
1114 for (v = 0; v < Points->n_points; v++) {
1115 double dist2, tmpdist2;
1116 double x, y, z;
1117
1118 dist2 = thresh2 + thresh2;
1119 x = Points->x[v];
1120 y = Points->y[v];
1121 z = Points->z[v];
1122
1123 /* Box */
1124 rect.boundary[0] = Points->x[v] - thresh;
1125 rect.boundary[3] = Points->x[v] + thresh;
1126 rect.boundary[1] = Points->y[v] - thresh;
1127 rect.boundary[4] = Points->y[v] + thresh;
1128 if (with_z) {
1129 rect.boundary[2] = Points->z[v] - thresh;
1130 rect.boundary[5] = Points->z[v] + thresh;
1131 }
1132
1134
1135 RTreeSearch(pnt_tree, &rect, add_item_box, (void *)List);
1136
1137 for (i = 0; i < List->n_values; i++) {
1138 double dx = List->box[i].E - Points->x[v];
1139 double dy = List->box[i].N - Points->y[v];
1140 double dz = 0;
1141
1142 if (with_z)
1143 dz = List->box[i].T - Points->z[v];
1144
1145 tmpdist2 = dx * dx + dy * dy + dz * dz;
1146
1147 if (tmpdist2 < dist2) {
1148 dist2 = tmpdist2;
1149
1150 x = List->box[i].E;
1151 y = List->box[i].N;
1152 z = List->box[i].T;
1153 }
1154 }
1155
1157 Points->x[v] = x;
1158 Points->y[v] = y;
1159 Points->z[v] = z;
1160
1161 changed = 1;
1162 if (nsnapped)
1163 (*nsnapped)++;
1164 }
1165 }
1166
1167 /* go through all vertices of the line to snap */
1168 /* find nearest reference segment */
1169 for (v = 0; v < Points->n_points; v++) {
1170 double dist2, tmpdist2;
1171 double x, y, z;
1172
1173 dist2 = thresh2 + thresh2;
1174 x = Points->x[v];
1175 y = Points->y[v];
1176 z = Points->z[v];
1177
1178 /* Box */
1179 rect.boundary[0] = Points->x[v] - thresh;
1180 rect.boundary[3] = Points->x[v] + thresh;
1181 rect.boundary[1] = Points->y[v] - thresh;
1182 rect.boundary[4] = Points->y[v] + thresh;
1183 if (with_z) {
1184 rect.boundary[2] = Points->z[v] - thresh;
1185 rect.boundary[5] = Points->z[v] + thresh;
1186 }
1187
1189
1190 RTreeSearch(seg_tree, &rect, add_item_box, (void *)List);
1191
1192 for (i = 0; i < List->n_values; i++) {
1193 double x1, y1, z1, x2, y2, z2;
1194 double tmpx, tmpy, tmpz;
1195 int status;
1196
1197 segment = List->id[i];
1198
1199 if (XSegs[segment] & X1W) {
1200 x1 = List->box[i].W;
1201 x2 = List->box[i].E;
1202 }
1203 else {
1204 x1 = List->box[i].E;
1205 x2 = List->box[i].W;
1206 }
1207 if (XSegs[segment] & Y1S) {
1208 y1 = List->box[i].S;
1209 y2 = List->box[i].N;
1210 }
1211 else {
1212 y1 = List->box[i].N;
1213 y2 = List->box[i].S;
1214 }
1215 if (XSegs[segment] & Z1B) {
1216 z1 = List->box[i].B;
1217 z2 = List->box[i].T;
1218 }
1219 else {
1220 z1 = List->box[i].T;
1221 z2 = List->box[i].B;
1222 }
1223
1224 /* Check the distance */
1226 Points->x[v], Points->y[v], Points->z[v], x1, y1, z1, x2, y2,
1227 z2, with_z, &tmpx, &tmpy, &tmpz, NULL, &status);
1228
1229 if (tmpdist2 < dist2 && status == 0) {
1230 dist2 = tmpdist2;
1231
1232 x = tmpx;
1233 y = tmpy;
1234 z = tmpz;
1235 }
1236 }
1237
1239 Points->x[v] = x;
1240 Points->y[v] = y;
1241 Points->z[v] = z;
1242
1243 changed = 1;
1244 if (nsnapped)
1245 (*nsnapped)++;
1246 }
1247 }
1248
1250 G_free(XSegs);
1251
1252 /* go through all segments of the line to snap */
1253 /* find nearest reference vertex, add this vertex */
1254 for (v = 0; v < Points->n_points - 1; v++) {
1255 double x1, x2, y1, y2, z1, z2;
1256 double xmin, xmax, ymin, ymax, zmin, zmax;
1257
1258 x1 = Points->x[v];
1259 x2 = Points->x[v + 1];
1260 y1 = Points->y[v];
1261 y2 = Points->y[v + 1];
1262 if (with_z) {
1263 z1 = Points->z[v];
1264 z2 = Points->z[v + 1];
1265 }
1266 else {
1267 z1 = z2 = 0;
1268 }
1269
1270 Vect_append_point(NPoints, Points->x[v], Points->y[v], Points->z[v]);
1271
1272 /* Box */
1273 if (x1 <= x2) {
1274 xmin = x1;
1275 xmax = x2;
1276 }
1277 else {
1278 xmin = x2;
1279 xmax = x1;
1280 }
1281 if (y1 <= y2) {
1282 ymin = y1;
1283 ymax = y2;
1284 }
1285 else {
1286 ymin = y2;
1287 ymax = y1;
1288 }
1289 if (z1 <= z2) {
1290 zmin = z1;
1291 zmax = z2;
1292 }
1293 else {
1294 zmin = z2;
1295 zmax = z1;
1296 }
1297
1298 rect.boundary[0] = xmin - thresh;
1299 rect.boundary[3] = xmax + thresh;
1300 rect.boundary[1] = ymin - thresh;
1301 rect.boundary[4] = ymax + thresh;
1302 rect.boundary[2] = zmin - thresh;
1303 rect.boundary[5] = zmax + thresh;
1304
1305 /* Find points */
1307 RTreeSearch(pnt_tree, &rect, add_item_box, (void *)List);
1308
1309 G_debug(3, " %d points in box", List->n_values);
1310
1311 /* Snap to vertex in threshold different from end points */
1312 nnew = 0;
1313 for (i = 0; i < List->n_values; i++) {
1314 double dist2, along;
1315 int status;
1316
1317 if (!with_z)
1318 List->box[i].T = 0;
1319
1320 if (Points->x[v] == List->box[i].E &&
1321 Points->y[v] == List->box[i].N &&
1322 Points->z[v] == List->box[i].T)
1323 continue; /* start point */
1324
1325 if (Points->x[v + 1] == List->box[i].E &&
1326 Points->y[v + 1] == List->box[i].N &&
1327 Points->z[v + 1] == List->box[i].T)
1328 continue; /* end point */
1329
1330 /* Check the distance */
1332 List->box[i].E, List->box[i].N, List->box[i].T, x1, y1, z1, x2,
1333 y2, z2, with_z, NULL, NULL, NULL, &along, &status);
1334
1335 if (dist2 <= thresh2 && status == 0) {
1336 G_debug(4, " anchor in thresh, along = %lf", along);
1337
1338 if (nnew == anew) {
1339 anew += 100;
1340 New = (NEW2 *)G_realloc(New, anew * sizeof(NEW2));
1341 }
1342 New[nnew].x = List->box[i].E;
1343 New[nnew].y = List->box[i].N;
1344 New[nnew].z = List->box[i].T;
1345 New[nnew].along = along;
1346 nnew++;
1347 }
1348 G_debug(3, "dist: %g, thresh: %g", dist2, thresh2);
1349 }
1350 G_debug(3, " nnew = %d", nnew);
1351 /* insert new vertices */
1352 if (nnew > 0) {
1353 /* sort by distance along the segment */
1354 qsort(New, nnew, sizeof(NEW2), sort_new2);
1355
1356 for (i = 0; i < nnew; i++) {
1357 Vect_append_point(NPoints, New[i].x, New[i].y, New[i].z);
1358 if (ncreated)
1359 (*ncreated)++;
1360 }
1361 changed = 1;
1362 }
1363 }
1364
1365 /* append end point */
1366 v = Points->n_points - 1;
1367 Vect_append_point(NPoints, Points->x[v], Points->y[v], Points->z[v]);
1368
1369 if (Points->n_points != NPoints->n_points) {
1370 Vect_line_prune(NPoints); /* remove duplicates */
1371 Vect_reset_line(Points);
1373 }
1374
1379 G_free(New);
1381 G_free(rect.boundary);
1382
1383 return changed;
1384}
#define X1W
Definition Vlib/snap.c:30
void Vect_snap_lines_list(struct Map_info *Map, const struct ilist *List_lines, double thresh, struct Map_info *Err)
Snap selected lines to existing vertex in threshold.
Definition Vlib/snap.c:172
#define Y1S
Definition Vlib/snap.c:31
#define Z1B
Definition Vlib/snap.c:32
void Vect_snap_lines(struct Map_info *Map, int type, double thresh, struct Map_info *Err)
Snap lines in vector map to existing vertex in threshold.
Definition Vlib/snap.c:909
int Vect_snap_line(struct Map_info *Map, struct ilist *reflist, struct line_pnts *Points, double thresh, int with_z, int *nsnapped, int *ncreated)
Snap a line to reference lines in Map with threshold.
Definition Vlib/snap.c:962
#define NULL
Definition ccmath.h:32
void G_percent(long, long, int)
Print percent complete messages.
Definition percent.c:61
void G_free(void *)
Free allocated memory.
Definition gis/alloc.c:147
#define G_realloc(p, n)
Definition defs/gis.h:141
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
#define G_malloc(n)
Definition defs/gis.h:139
void void G_verbose_message(const char *,...) __attribute__((format(printf
char * G_tempfile(void)
Returns a temporary file name.
Definition tempfile.c:62
void void void G_important_message(const char *,...) __attribute__((format(printf
void G_ilist_add(struct ilist *, int)
Add item to ilist.
Definition ilist.c:78
int G_debug(int, const char *,...) __attribute__((format(printf
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
Definition line.c:77
off_t Vect_rewrite_line(struct Map_info *, off_t, int, const struct line_pnts *, const struct line_cats *)
Rewrites existing feature (topological level required)
int Vect_reset_boxlist(struct boxlist *)
Reset boxlist structure.
plus_t Vect_get_num_lines(struct Map_info *)
Fetch number of features (points, lines, boundaries, centroids) in vector map.
Definition level_two.c:75
struct boxlist * Vect_new_boxlist(int)
Creates and initializes a struct boxlist.
void Vect_destroy_boxlist(struct boxlist *)
Frees all memory associated with a struct boxlist, including the struct itself.
void Vect_destroy_list(struct ilist *)
Frees all memory associated with a struct ilist, including the struct itself.
void Vect_destroy_cats_struct(struct line_cats *)
Frees all memory associated with line_cats structure, including the struct itself.
int Vect_read_line(struct Map_info *, struct line_pnts *, struct line_cats *, int)
Read vector feature (topological level required)
int Vect_line_alive(struct Map_info *, int)
Check if feature is alive or dead (topological level required)
int Vect_delete_line(struct Map_info *, off_t)
Delete existing feature (topological level required)
struct line_cats * Vect_new_cats_struct(void)
Creates and initializes line_cats structure.
struct ilist * Vect_new_list(void)
Creates and initializes a struct ilist.
off_t Vect_write_line(struct Map_info *, int, const struct line_pnts *, const struct line_cats *)
Writes a new feature.
void Vect_reset_line(struct line_pnts *)
Reset line.
Definition line.c:129
int Vect_line_prune(struct line_pnts *)
Remove duplicate points, i.e. zero length segments.
Definition line.c:279
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition line.c:45
int Vect_reset_list(struct ilist *)
Reset ilist structure.
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition line.c:148
int Vect_append_points(struct line_pnts *, const struct line_pnts *, int)
Appends points to the end of a line.
Definition line.c:335
#define GV_LINES
#define GV_FORWARD
Line direction indicator forward/backward.
double dig_distance2_point_to_line(double, double, double, double, double, double, double, double, double, int, double *, double *, double *, double *, int *)
int dig_boxlist_add(struct boxlist *, int, const struct bound_box *)
Header file for msvc/fcntl.c.
#define open
Definition fcntl.h:33
#define UNUSED
A macro for an attribute, if attached to a variable, indicating that the variable is not used.
Definition gis.h:46
#define _(str)
Definition glocale.h:10
int kdtree_rnn(struct kdtree *t, double *c, int **puid, int *skip)
Definition kdtree.c:751
struct kdtree * kdtree_create(char ndims, int *btol)
Definition kdtree.c:111
int kdtree_insert(struct kdtree *t, double *c, int uid, int dc)
Definition kdtree.c:179
int kdtree_knn(struct kdtree *t, double *c, int *uid, double *d, int k, int *skip)
Definition kdtree.c:512
void kdtree_destroy(struct kdtree *t)
Definition kdtree.c:167
int kdtree_dnn(struct kdtree *t, double *c, int **puid, double **pd, double maxdist, int *skip)
Definition kdtree.c:636
Vector map info.
RectReal * boundary
Definition rtree.h:55
Definition rtree.h:123
Bounding box.
Definition dig_structs.h:62
double W
West.
Definition dig_structs.h:78
List of bounding boxes with id.
List of integers.
Definition gis.h:715
k-d tree
Definition kdtree.h:80
Feature category info.
Feature geometry info - coordinates.
double * y
Array of Y coordinates.
double * x
Array of X coordinates.
int n_points
Number of points.
double * z
Array of Z coordinates.
Definition manage.h:4
#define close
Definition unistd.h:8
void RTreeSetOverflow(struct RTree *t, char overflow)
Enable/disable R*-tree forced reinsertion (overflow)
struct RTree * RTreeCreateTree(int fd, off_t rootpos, int ndims)
Create new empty R*-Tree.
int RTreeInsertRect(struct RTree_Rect *r, int tid, struct RTree *t)
Insert an item into a R*-Tree.
void RTreeDestroyTree(struct RTree *t)
Destroy an R*-Tree.
int RTreeSearch(struct RTree *t, struct RTree_Rect *r, SearchHitCallback *shcb, void *cbarg)
Search an R*-Tree.
#define x