24 #include <grass/kdtree.h> 47 static int sort_new(
const void *pa,
const void *pb)
52 return (p1->along < p2->along ? -1 : (p1->along > p2->along));
65 double x, y, z, along;
69 static int sort_new2(
const void *pa,
const void *pb)
71 NEW2 *p1 = (NEW2 *) pa;
72 NEW2 *p2 = (NEW2 *) pb;
74 return (p1->along < p2->along ? -1 : (p1->along > p2->along));
78 static int find_item(
int id,
const struct RTree_Rect *rect,
void *
list)
85 static int add_item(
int id,
const struct RTree_Rect *rect,
void *
list)
92 static int find_item_box(
int id,
const struct RTree_Rect *rect,
void *
list)
109 static int add_item_box(
int id,
const struct RTree_Rect *rect,
void *
list)
126 Vect_snap_lines_list_rtree(
struct Map_info *,
const struct ilist *,
130 Vect_snap_lines_list_kdtree(
struct Map_info *,
const struct ilist *,
171 double thresh,
struct Map_info *Err)
173 if (
getenv(
"GRASS_VECTOR_LOWMEM"))
174 Vect_snap_lines_list_rtree(Map, List_lines, thresh, Err);
176 Vect_snap_lines_list_kdtree(Map, List_lines, thresh, Err);
180 Vect_snap_lines_list_kdtree(
struct Map_info *Map,
const struct ilist *List_lines,
181 double thresh,
struct Map_info *Err)
185 int line, ltype, line_idx;
189 int nanchors, ntosnap;
190 int nsnapped, ncreated;
191 int apoints, npoints, nvertices;
202 int *kduid, kd_found;
215 thresh2 = thresh * thresh;
224 for (line_idx = 0; line_idx < List_lines->
n_values; line_idx++) {
229 line = List_lines->
value[line_idx];
237 for (v = 0; v < Points->
n_points; v++) {
239 G_debug(3,
" vertex v = %d", v);
248 if ((point - 1) == apoints) {
252 (apoints + 1) *
sizeof(XPNT));
254 XPnts[point].x = Points->
x[v];
255 XPnts[point].y = Points->
y[v];
256 XPnts[point].anchor = -1;
270 nanchors = ntosnap = 0;
271 for (point = 1; point <= npoints; point++) {
276 G_debug(3,
" point = %d", point);
278 if (XPnts[point].anchor >= 0)
281 XPnts[point].anchor = 0;
285 c[0] = XPnts[point].x;
286 c[1] = XPnts[point].y;
289 kd_found =
kdtree_dnn(KDTree, c, &kduid, &kdd, thresh, &point);
290 G_debug(4,
" %d points in threshold box", kd_found);
292 for (i = 0; i < kd_found; i++) {
294 double dx, dy, dist2;
300 dx = XPnts[pointb].x - XPnts[point].x;
301 dy = XPnts[pointb].y - XPnts[point].y;
302 dist2 = dx * dx + dy * dy;
308 if (XPnts[pointb].anchor == -1) {
309 XPnts[pointb].anchor = point;
312 else if (XPnts[pointb].anchor > 0) {
315 dx = XPnts[XPnts[pointb].anchor].x - XPnts[pointb].x;
316 dy = XPnts[XPnts[pointb].anchor].y - XPnts[pointb].y;
317 dist2_a = dx * dx + dy * dy;
320 if (dist2 < dist2_a) {
321 XPnts[pointb].anchor = point;
335 nsnapped = ncreated = 0;
339 for (line_idx = 0; line_idx < List_lines->
n_values; line_idx++) {
340 int v, spoint, anchor;
346 line = List_lines->
value[line_idx];
356 Index = (
int *)
G_realloc(Index, aindex *
sizeof(
int));
360 G_debug(3,
"Snap all vertices");
361 for (v = 0; v < Points->
n_points; v++) {
374 anchor = XPnts[spoint].anchor;
377 Points->
x[v] = XPnts[anchor].x;
378 Points->
y[v] = XPnts[anchor].y;
392 G_debug(3,
"Snap all segments");
393 for (v = 0; v < Points->
n_points - 1; v++) {
395 double x1, x2, y1, y2, xmin, xmax, ymin, ymax;
398 G_debug(3,
" segment = %d end anchors : %d %d", v, Index[v],
402 x2 = Points->
x[v + 1];
404 y2 = Points->
y[v + 1];
429 G_debug(3,
" search anchors for segment %g,%g to %g,%g", x1, y1, x2, y2);
434 rc[0] = xmin - thresh * 2;
435 rc[1] = ymin - thresh * 2;
436 rc[2] = xmax + thresh * 2;
437 rc[3] = ymax + thresh * 2;
441 G_debug(3,
" %d points in box", kd_found);
445 for (i = 0; i < kd_found; i++) {
450 G_debug(4,
" spoint = %d anchor = %d", spoint,
451 XPnts[spoint].anchor);
453 if (spoint == Index[v] || spoint == Index[v + 1])
455 if (XPnts[spoint].anchor > 0)
461 XPnts[spoint].y, 0, x1, y1, 0,
463 NULL, &along, &status);
465 G_debug(4,
" distance = %lf", sqrt(dist2));
467 if (status == 0 && dist2 <= thresh2) {
468 G_debug(4,
" anchor in thresh, along = %lf", along);
472 New = (NEW *)
G_realloc(New, anew *
sizeof(NEW));
474 New[nnew].anchor = spoint;
475 New[nnew].along = along;
482 G_debug(3,
" nnew = %d", nnew);
486 qsort(New,
sizeof(
char) * nnew,
sizeof(NEW), sort_new);
488 for (i = 0; i < nnew; i++) {
489 anchor = New[i].anchor;
531 Vect_snap_lines_list_rtree(
struct Map_info *Map,
const struct ilist *List_lines,
532 double thresh,
struct Map_info *Err)
536 int line, ltype, line_idx;
540 int nanchors, ntosnap;
541 int nsnapped, ncreated;
542 int apoints, npoints, nvertices;
553 static int rect_init = 0;
567 if (
getenv(
"GRASS_VECTOR_LOWMEM")) {
570 rtreefd = open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
575 thresh2 = thresh * thresh;
584 for (line_idx = 0; line_idx < List_lines->
n_values; line_idx++) {
589 line = List_lines->
value[line_idx];
597 for (v = 0; v < Points->
n_points; v++) {
598 G_debug(3,
" vertex v = %d", v);
617 if ((point - 1) == apoints) {
621 (apoints + 1) *
sizeof(XPNT));
623 XPnts[point].x = Points->
x[v];
624 XPnts[point].y = Points->
y[v];
625 XPnts[point].anchor = -1;
639 nanchors = ntosnap = 0;
640 for (point = 1; point <= npoints; point++) {
645 G_debug(3,
" point = %d", point);
647 if (XPnts[point].anchor >= 0)
650 XPnts[point].anchor = 0;
654 rect.
boundary[0] = XPnts[point].x - thresh;
655 rect.
boundary[3] = XPnts[point].x + thresh;
656 rect.
boundary[1] = XPnts[point].y - thresh;
657 rect.
boundary[4] = XPnts[point].y + thresh;
665 for (i = 0; i < List->
n_values; i++) {
667 double dx, dy, dist2;
669 pointb = List->
value[i];
673 dx = XPnts[pointb].x - XPnts[point].x;
674 dy = XPnts[pointb].y - XPnts[point].y;
675 dist2 = dx * dx + dy * dy;
681 if (XPnts[pointb].anchor == -1) {
682 XPnts[pointb].anchor = point;
685 else if (XPnts[pointb].anchor > 0) {
688 dx = XPnts[XPnts[pointb].anchor].x - XPnts[pointb].x;
689 dy = XPnts[XPnts[pointb].anchor].y - XPnts[pointb].y;
690 dist2_a = dx * dx + dy * dy;
693 if (dist2 < dist2_a) {
694 XPnts[pointb].anchor = point;
704 nsnapped = ncreated = 0;
708 for (line_idx = 0; line_idx < List_lines->
n_values; line_idx++) {
709 int v, spoint, anchor;
714 line = List_lines->
value[line_idx];
724 Index = (
int *)
G_realloc(Index, aindex *
sizeof(
int));
728 for (v = 0; v < Points->
n_points; v++) {
742 spoint = List->
value[0];
743 anchor = XPnts[spoint].anchor;
746 Points->
x[v] = XPnts[anchor].x;
747 Points->
y[v] = XPnts[anchor].y;
761 for (v = 0; v < Points->
n_points - 1; v++) {
763 double x1, x2, y1, y2, xmin, xmax, ymin, ymax;
765 G_debug(3,
" segment = %d end anchors : %d %d", v, Index[v],
769 x2 = Points->
x[v + 1];
771 y2 = Points->
y[v + 1];
809 for (i = 0; i < List->
n_values; i++) {
813 spoint = List->
value[i];
814 G_debug(4,
" spoint = %d anchor = %d", spoint,
815 XPnts[spoint].anchor);
817 if (spoint == Index[v] || spoint == Index[v + 1])
819 if (XPnts[spoint].anchor > 0)
825 XPnts[spoint].y, 0, x1, y1, 0,
827 NULL, &along, &status);
829 G_debug(4,
" distance = %lf", sqrt(dist2));
831 if (status == 0 && dist2 <= thresh2) {
832 G_debug(4,
" anchor in thresh, along = %lf", along);
836 New = (NEW *)
G_realloc(New, anew *
sizeof(NEW));
838 New[nnew].anchor = spoint;
839 New[nnew].along = along;
843 G_debug(3,
" nnew = %d", nnew);
847 qsort(New,
sizeof(
char) * nnew,
sizeof(NEW), sort_new);
849 for (i = 0; i < nnew; i++) {
850 anchor = New[i].anchor;
910 int line, nlines, ltype;
918 for (line = 1; line <= nlines; line++) {
962 struct line_pnts *Points,
double thresh,
int with_z,
963 int *nsnapped,
int *ncreated)
967 int i, v, line, nlines;
980 struct RTree *pnt_tree,
1011 with_z = (with_z != 0);
1017 thresh2 = thresh * thresh;
1019 point = segment = 1;
1026 for (i = 0; i < nlines; i++) {
1028 line = reflist->
value[i];
1030 G_debug(3,
"line = %d", line);
1037 for (v = 0; v < LPoints->
n_points; v++) {
1038 G_debug(3,
" vertex v = %d", v);
1053 RTreeSearch(pnt_tree, &rect, find_item_box, (
void *)List);
1054 G_debug(3,
"List : nvalues = %d", List->n_values);
1056 if (List->n_values == 0) {
1069 if (LPoints->
x[v - 1] < LPoints->
x[v]) {
1078 if (LPoints->
y[v - 1] < LPoints->
y[v]) {
1087 if (LPoints->
z[v - 1] < LPoints->
z[v]) {
1101 if ((segment - 1) == asegments) {
1105 (asegments + 1) *
sizeof(char));
1107 XSegs[segment] = sides;
1115 for (v = 0; v < Points->
n_points; v++) {
1116 double dist2, tmpdist2;
1119 dist2 = thresh2 + thresh2;
1125 rect.
boundary[0] = Points->
x[v] - thresh;
1126 rect.
boundary[3] = Points->
x[v] + thresh;
1127 rect.
boundary[1] = Points->
y[v] - thresh;
1128 rect.
boundary[4] = Points->
y[v] + thresh;
1130 rect.
boundary[2] = Points->
z[v] - thresh;
1131 rect.
boundary[5] = Points->
z[v] + thresh;
1136 RTreeSearch(pnt_tree, &rect, add_item_box, (
void *)List);
1138 for (i = 0; i < List->n_values; i++) {
1139 double dx = List->box[i].E - Points->
x[v];
1140 double dy = List->box[i].N - Points->
y[v];
1144 dz = List->box[i].T - Points->
z[v];
1146 tmpdist2 = dx * dx + dy * dy + dz * dz;
1148 if (tmpdist2 < dist2) {
1157 if (dist2 <= thresh2 && dist2 > 0) {
1170 for (v = 0; v < Points->
n_points; v++) {
1171 double dist2, tmpdist2;
1174 dist2 = thresh2 + thresh2;
1180 rect.
boundary[0] = Points->
x[v] - thresh;
1181 rect.
boundary[3] = Points->
x[v] + thresh;
1182 rect.
boundary[1] = Points->
y[v] - thresh;
1183 rect.
boundary[4] = Points->
y[v] + thresh;
1185 rect.
boundary[2] = Points->
z[v] - thresh;
1186 rect.
boundary[5] = Points->
z[v] + thresh;
1191 RTreeSearch(seg_tree, &rect, add_item_box, (
void *)List);
1193 for (i = 0; i < List->n_values; i++) {
1194 double x1, y1, z1, x2, y2, z2;
1195 double tmpx, tmpy, tmpz;
1198 segment = List->id[i];
1200 if (XSegs[segment] &
X1W) {
1201 x1 = List->box[i].W;
1202 x2 = List->box[i].E;
1205 x1 = List->box[i].E;
1206 x2 = List->box[i].W;
1208 if (XSegs[segment] &
Y1S) {
1209 y1 = List->box[i].S;
1210 y2 = List->box[i].N;
1213 y1 = List->box[i].N;
1214 y2 = List->box[i].S;
1216 if (XSegs[segment] &
Z1B) {
1217 z1 = List->box[i].B;
1218 z2 = List->box[i].T;
1221 z1 = List->box[i].T;
1222 z2 = List->box[i].B;
1230 x1, y1, z1, x2, y2, z2,
1231 with_z, &tmpx, &tmpy, &tmpz,
1234 if (tmpdist2 < dist2 && status == 0) {
1243 if (dist2 <= thresh2 && dist2 > 0) {
1259 for (v = 0; v < Points->
n_points - 1; v++) {
1260 double x1, x2, y1, y2, z1, z2;
1261 double xmin, xmax, ymin, ymax, zmin, zmax;
1264 x2 = Points->
x[v + 1];
1266 y2 = Points->
y[v + 1];
1269 z2 = Points->
z[v + 1];
1313 RTreeSearch(pnt_tree, &rect, add_item_box, (
void *)List);
1315 G_debug(3,
" %d points in box", List->n_values);
1319 for (i = 0; i < List->n_values; i++) {
1320 double dist2, along;
1326 if (Points->
x[v] == List->box[i].E &&
1327 Points->
y[v] == List->box[i].N &&
1328 Points->
z[v] == List->box[i].T)
1331 if (Points->
x[v + 1] == List->box[i].E &&
1332 Points->
y[v + 1] == List->box[i].N &&
1333 Points->
z[v + 1] == List->box[i].T)
1340 List->box[i].T, x1, y1, z1,
1342 NULL, &along, &status);
1344 if (dist2 <= thresh2 && status == 0) {
1345 G_debug(4,
" anchor in thresh, along = %lf", along);
1349 New = (NEW2 *)
G_realloc(New, anew *
sizeof(NEW2));
1351 New[nnew].x = List->box[i].E;
1352 New[nnew].y = List->box[i].N;
1353 New[nnew].z = List->box[i].T;
1354 New[nnew].along = along;
1357 G_debug(3,
"dist: %g, thresh: %g", dist2, thresh2);
1359 G_debug(3,
" nnew = %d", nnew);
1363 qsort(New, nnew,
sizeof(NEW2), sort_new2);
1365 for (i = 0; i < nnew; i++) {
void Vect_destroy_boxlist(struct boxlist *)
Frees all memory associated with a struct boxlist, including the struct itself.
int Vect_append_points(struct line_pnts *, const struct line_pnts *, int)
Appends points to the end of a line.
int RTreeInsertRect(struct RTree_Rect *r, int tid, struct RTree *t)
Insert an item into a R*-Tree.
#define GV_FORWARD
Line direction indicator forward/backward.
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void void void G_important_message(const char *,...) __attribute__((format(printf
int Vect_reset_list(struct ilist *)
Reset ilist structure.
void RTreeDestroyTree(struct RTree *t)
Destroy an R*-Tree.
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.
plus_t Vect_get_num_lines(const struct Map_info *)
Fetch number of features (points, lines, boundaries, centroids) in vector map.
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 RTreeSearch(struct RTree *t, struct RTree_Rect *r, SearchHitCallback *shcb, void *cbarg)
Search an R*-Tree.
int n_points
Number of points.
int kdtree_rnn(struct kdtree *t, double *c, int **puid, int *skip)
int n_values
Number of values in the list.
void G_free(void *)
Free allocated memory.
int kdtree_dnn(struct kdtree *t, double *c, int **puid, double **pd, double maxdist, int *skip)
int Vect_line_prune(struct line_pnts *)
Remove duplicate points, i.e. zero length segments.
struct kdtree * kdtree_create(char ndims, int *btol)
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.
int dig_boxlist_add(struct boxlist *, int, const struct bound_box *)
double * x
Array of X coordinates.
void G_ilist_add(struct ilist *, int)
Add item to ilist.
Feature geometry info - coordinates.
void Vect_destroy_list(struct ilist *)
Frees all memory associated with a struct ilist, including the struct itself.
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
void Vect_destroy_cats_struct(struct line_cats *)
Frees all memory associated with line_cats structure, including the struct itself.
char * G_tempfile(void)
Returns a temporary file name.
off_t Vect_write_line(struct Map_info *, int, const struct line_pnts *, const struct line_cats *)
Writes a new feature.
double dig_distance2_point_to_line(double, double, double, double, double, double, double, double, double, int, double *, double *, double *, double *, int *)
void RTreeSetOverflow(struct RTree *t, char overflow)
Enable/disable R*-tree forced reinsertion (overflow)
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
void kdtree_destroy(struct kdtree *t)
int kdtree_knn(struct kdtree *t, double *c, int *uid, double *d, int k, int *skip)
int kdtree_insert(struct kdtree *t, double *c, int uid, int dc)
int Vect_reset_boxlist(struct boxlist *)
Reset boxlist structure.
void G_percent(long, long, int)
Print percent complete messages.
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.
double * y
Array of Y coordinates.
struct ilist * Vect_new_list(void)
Creates and initializes a struct ilist.
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.
List of bounding boxes with id.
int Vect_line_alive(const struct Map_info *, int)
Check if feature is alive or dead (topological level required)
double * z
Array of Z coordinates.
int * value
Array of values.
struct boxlist * Vect_new_boxlist(int)
Creates and initializes a struct boxlist.
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
void void G_verbose_message(const char *,...) __attribute__((format(printf
int Vect_read_line(const struct Map_info *, struct line_pnts *, struct line_cats *, int)
Read vector feature (topological level required)
int G_debug(int, const char *,...) __attribute__((format(printf
struct RTree * RTreeCreateTree(int fd, off_t rootpos, int ndims)
Create new empty R*-Tree.
void Vect_reset_line(struct line_pnts *)
Reset line.