27 #include <grass/kdtree.h>
49 static int sort_new(
const void *pa,
const void *pb)
54 return (p1->along < p2->along ? -1 : (p1->along > p2->along));
66 double x, y, z, along;
70 static int sort_new2(
const void *pa,
const void *pb)
72 NEW2 *p1 = (NEW2 *)pa;
73 NEW2 *p2 = (NEW2 *)pb;
75 return (p1->along < p2->along ? -1 : (p1->along > p2->along));
95 static int find_item_box(
int id,
const struct RTree_Rect *rect,
void *
list)
113 static int add_item_box(
int id,
const struct RTree_Rect *rect,
void *
list)
129 static void Vect_snap_lines_list_rtree(
struct Map_info *,
const struct ilist *,
132 static void Vect_snap_lines_list_kdtree(
struct Map_info *,
const struct ilist *,
173 double thresh,
struct Map_info *Err)
175 if (getenv(
"GRASS_VECTOR_LOWMEM"))
176 Vect_snap_lines_list_rtree(Map, List_lines, thresh, Err);
178 Vect_snap_lines_list_kdtree(Map, List_lines, thresh, Err);
181 static void Vect_snap_lines_list_kdtree(
struct Map_info *Map,
182 const struct ilist *List_lines,
183 double thresh,
struct Map_info *Err)
187 int line, ltype, line_idx;
191 int nsnapped, ncreated;
193 int apoints, npoints;
204 int *kduid, kd_found;
216 thresh2 = thresh * thresh;
225 for (line_idx = 0; line_idx < List_lines->
n_values; line_idx++) {
230 line = List_lines->
value[line_idx];
238 for (v = 0; v < Points->
n_points; v++) {
240 G_debug(3,
" vertex v = %d", v);
248 if ((point - 1) == apoints) {
251 (XPNT *)
G_realloc(XPnts, (apoints + 1) *
sizeof(XPNT));
253 XPnts[point].x = Points->
x[v];
254 XPnts[point].y = Points->
y[v];
255 XPnts[point].anchor = -1;
269 for (point = 1; point <= npoints; point++) {
274 G_debug(3,
" point = %d", point);
276 if (XPnts[point].anchor >= 0)
279 XPnts[point].anchor = 0;
282 c[0] = XPnts[point].x;
283 c[1] = XPnts[point].y;
286 kd_found =
kdtree_dnn(KDTree, c, &kduid, &kdd, thresh, &point);
287 G_debug(4,
" %d points in threshold box", kd_found);
289 for (i = 0; i < kd_found; i++) {
291 double dx, dy, dist2;
297 dx = XPnts[pointb].x - XPnts[point].x;
298 dy = XPnts[pointb].y - XPnts[point].y;
299 dist2 = dx * dx + dy * dy;
305 if (XPnts[pointb].anchor == -1) {
306 XPnts[pointb].anchor = point;
308 else if (XPnts[pointb].anchor >
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;
317 if (dist2 < dist2_a) {
318 XPnts[pointb].anchor = point;
333 nsnapped = ncreated = 0;
337 for (line_idx = 0; line_idx < List_lines->
n_values; line_idx++) {
338 int v, spoint, anchor;
344 line = List_lines->
value[line_idx];
354 Index = (
int *)
G_realloc(Index, aindex *
sizeof(
int));
358 G_debug(3,
"Snap all vertices");
359 for (v = 0; v < Points->
n_points; v++) {
372 anchor = XPnts[spoint].anchor;
375 Points->
x[v] = XPnts[anchor].x;
376 Points->
y[v] = XPnts[anchor].y;
390 G_debug(3,
"Snap all segments");
391 for (v = 0; v < Points->
n_points - 1; v++) {
393 double x1, x2, y1, y2, xmin, xmax, ymin, ymax;
396 G_debug(3,
" segment = %d end anchors : %d %d", v, Index[v],
400 x2 = Points->
x[v + 1];
402 y2 = Points->
y[v + 1];
427 G_debug(3,
" search anchors for segment %g,%g to %g,%g", x1, y1,
433 rc[0] = xmin - thresh * 2;
434 rc[1] = ymin - thresh * 2;
435 rc[2] = xmax + thresh * 2;
436 rc[3] = ymax + thresh * 2;
440 G_debug(3,
" %d points in box", kd_found);
444 for (i = 0; i < kd_found; i++) {
449 G_debug(4,
" spoint = %d anchor = %d", spoint,
450 XPnts[spoint].anchor);
452 if (spoint == Index[v] || spoint == Index[v + 1])
454 if (XPnts[spoint].anchor > 0)
459 XPnts[spoint].
x, XPnts[spoint].y, 0, x1, y1, 0, x2, y2, 0,
462 G_debug(4,
" distance = %lf", sqrt(dist2));
464 if (status == 0 && dist2 <= thresh2) {
465 G_debug(4,
" anchor in thresh, along = %lf", along);
469 New = (NEW *)
G_realloc(New, anew *
sizeof(NEW));
471 New[nnew].anchor = spoint;
472 New[nnew].along = along;
479 G_debug(3,
" nnew = %d", nnew);
483 qsort(New,
sizeof(
char) * nnew,
sizeof(NEW), sort_new);
485 for (i = 0; i < nnew; i++) {
486 anchor = New[i].anchor;
529 static void Vect_snap_lines_list_rtree(
struct Map_info *Map,
530 const struct ilist *List_lines,
531 double thresh,
struct Map_info *Err)
535 int line, ltype, line_idx;
539 int nsnapped, ncreated;
541 int apoints, npoints;
552 static int rect_init = 0;
566 if (getenv(
"GRASS_VECTOR_LOWMEM")) {
569 rtreefd =
open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
572 filename, strerror(errno));
574 if (remove(filename) != 0) {
575 G_warning(
_(
"Unable to remove temporary file <%s>: %s"), filename,
582 thresh2 = thresh * thresh;
591 for (line_idx = 0; line_idx < List_lines->
n_values; line_idx++) {
596 line = List_lines->
value[line_idx];
604 for (v = 0; v < Points->
n_points; v++) {
605 G_debug(3,
" vertex v = %d", v);
623 if ((point - 1) == apoints) {
626 (XPNT *)
G_realloc(XPnts, (apoints + 1) *
sizeof(XPNT));
628 XPnts[point].x = Points->
x[v];
629 XPnts[point].y = Points->
y[v];
630 XPnts[point].anchor = -1;
644 for (point = 1; point <= npoints; point++) {
649 G_debug(3,
" point = %d", point);
651 if (XPnts[point].anchor >= 0)
654 XPnts[point].anchor = 0;
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;
668 for (i = 0; i < List->
n_values; i++) {
670 double dx, dy, dist2;
672 pointb = List->
value[i];
676 dx = XPnts[pointb].x - XPnts[point].x;
677 dy = XPnts[pointb].y - XPnts[point].y;
678 dist2 = dx * dx + dy * dy;
684 if (XPnts[pointb].anchor == -1) {
685 XPnts[pointb].anchor = point;
687 else if (XPnts[pointb].anchor >
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;
696 if (dist2 < dist2_a) {
697 XPnts[pointb].anchor = point;
708 nsnapped = ncreated = 0;
712 for (line_idx = 0; line_idx < List_lines->
n_values; line_idx++) {
713 int v, spoint, anchor;
718 line = List_lines->
value[line_idx];
728 Index = (
int *)
G_realloc(Index, aindex *
sizeof(
int));
732 for (v = 0; v < Points->
n_points; v++) {
746 spoint = List->
value[0];
747 anchor = XPnts[spoint].anchor;
750 Points->
x[v] = XPnts[anchor].x;
751 Points->
y[v] = XPnts[anchor].y;
765 for (v = 0; v < Points->
n_points - 1; v++) {
767 double x1, x2, y1, y2, xmin, xmax, ymin, ymax;
769 G_debug(3,
" segment = %d end anchors : %d %d", v, Index[v],
773 x2 = Points->
x[v + 1];
775 y2 = Points->
y[v + 1];
813 for (i = 0; i < List->
n_values; i++) {
817 spoint = List->
value[i];
818 G_debug(4,
" spoint = %d anchor = %d", spoint,
819 XPnts[spoint].anchor);
821 if (spoint == Index[v] || spoint == Index[v + 1])
823 if (XPnts[spoint].anchor > 0)
828 XPnts[spoint].
x, XPnts[spoint].y, 0, x1, y1, 0, x2, y2, 0,
831 G_debug(4,
" distance = %lf", sqrt(dist2));
833 if (status == 0 && dist2 <= thresh2) {
834 G_debug(4,
" anchor in thresh, along = %lf", along);
838 New = (NEW *)
G_realloc(New, anew *
sizeof(NEW));
840 New[nnew].anchor = spoint;
841 New[nnew].along = along;
845 G_debug(3,
" nnew = %d", nnew);
849 qsort(New,
sizeof(
char) * nnew,
sizeof(NEW), sort_new);
851 for (i = 0; i < nnew; i++) {
852 anchor = New[i].anchor;
912 int line, nlines, ltype;
920 for (line = 1; line <= nlines; line++) {
963 struct line_pnts *Points,
double thresh,
int with_z,
964 int *nsnapped,
int *ncreated)
968 int i, v, line, nlines;
980 struct RTree *pnt_tree,
1013 with_z = (with_z != 0);
1019 thresh2 = thresh * thresh;
1021 point = segment = 1;
1027 for (i = 0; i < nlines; i++) {
1029 line = reflist->
value[i];
1031 G_debug(3,
"line = %d", line);
1038 for (v = 0; v < LPoints->
n_points; v++) {
1039 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) {
1104 (asegments + 1) *
sizeof(char));
1106 XSegs[segment] = sides;
1114 for (v = 0; v < Points->
n_points; v++) {
1115 double dist2, tmpdist2;
1118 dist2 = thresh2 + thresh2;
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;
1129 rect.
boundary[2] = Points->
z[v] - thresh;
1130 rect.
boundary[5] = Points->
z[v] + thresh;
1135 RTreeSearch(pnt_tree, &rect, add_item_box, (
void *)List);
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];
1143 dz = List->box[i].T - Points->
z[v];
1145 tmpdist2 = dx * dx + dy * dy + dz * dz;
1147 if (tmpdist2 < dist2) {
1156 if (dist2 <= thresh2 && dist2 > 0) {
1169 for (v = 0; v < Points->
n_points; v++) {
1170 double dist2, tmpdist2;
1173 dist2 = thresh2 + thresh2;
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;
1184 rect.
boundary[2] = Points->
z[v] - thresh;
1185 rect.
boundary[5] = Points->
z[v] + thresh;
1190 RTreeSearch(seg_tree, &rect, add_item_box, (
void *)List);
1192 for (i = 0; i < List->n_values; i++) {
1193 double x1, y1, z1, x2, y2, z2;
1194 double tmpx, tmpy, tmpz;
1197 segment = List->id[i];
1199 if (XSegs[segment] &
X1W) {
1200 x1 = List->box[i].W;
1201 x2 = List->box[i].E;
1204 x1 = List->box[i].E;
1205 x2 = List->box[i].W;
1207 if (XSegs[segment] &
Y1S) {
1208 y1 = List->box[i].S;
1209 y2 = List->box[i].N;
1212 y1 = List->box[i].N;
1213 y2 = List->box[i].S;
1215 if (XSegs[segment] &
Z1B) {
1216 z1 = List->box[i].B;
1217 z2 = List->box[i].T;
1220 z1 = List->box[i].T;
1221 z2 = List->box[i].B;
1226 Points->
x[v], Points->
y[v], Points->
z[v], x1, y1, z1, x2, y2,
1227 z2, with_z, &tmpx, &tmpy, &tmpz,
NULL, &status);
1229 if (tmpdist2 < dist2 && status == 0) {
1238 if (dist2 <= thresh2 && dist2 > 0) {
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;
1259 x2 = Points->
x[v + 1];
1261 y2 = Points->
y[v + 1];
1264 z2 = Points->
z[v + 1];
1307 RTreeSearch(pnt_tree, &rect, add_item_box, (
void *)List);
1309 G_debug(3,
" %d points in box", List->n_values);
1313 for (i = 0; i < List->n_values; i++) {
1314 double dist2, along;
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)
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)
1332 List->box[i].E, List->box[i].N, List->box[i].T, x1, y1, z1, x2,
1335 if (dist2 <= thresh2 && status == 0) {
1336 G_debug(4,
" anchor in thresh, along = %lf", along);
1340 New = (NEW2 *)
G_realloc(New, anew *
sizeof(NEW2));
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;
1348 G_debug(3,
"dist: %g, thresh: %g", dist2, thresh2);
1350 G_debug(3,
" nnew = %d", nnew);
1354 qsort(New, nnew,
sizeof(NEW2), sort_new2);
1356 for (i = 0; i < nnew; i++) {
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.
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 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.
void G_percent(long, long, int)
Print percent complete messages.
void G_free(void *)
Free allocated memory.
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
void void G_verbose_message(const char *,...) __attribute__((format(printf
char * G_tempfile(void)
Returns a temporary file name.
void void void G_important_message(const char *,...) __attribute__((format(printf
void G_ilist_add(struct ilist *, int)
Add item to ilist.
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.
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.
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)
struct ilist * Vect_new_list(void)
Creates and initializes a struct ilist.
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)
off_t Vect_write_line(struct Map_info *, int, const struct line_pnts *, const struct line_cats *)
Writes a new feature.
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
void Vect_reset_line(struct line_pnts *)
Reset line.
int Vect_line_prune(struct line_pnts *)
Remove duplicate points, i.e. zero length segments.
struct line_cats * Vect_new_cats_struct(void)
Creates and initializes line_cats structure.
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.
int Vect_append_points(struct line_pnts *, const struct line_pnts *, int)
Appends points to the end of a line.
#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 UNUSED
A macro for an attribute, if attached to a variable, indicating that the variable is not used.
int kdtree_rnn(struct kdtree *t, double *c, int **puid, int *skip)
struct kdtree * kdtree_create(char ndims, int *btol)
int kdtree_insert(struct kdtree *t, double *c, int uid, int dc)
int kdtree_knn(struct kdtree *t, double *c, int *uid, double *d, int k, int *skip)
void kdtree_destroy(struct kdtree *t)
int kdtree_dnn(struct kdtree *t, double *c, int **puid, double **pd, double maxdist, int *skip)
List of bounding boxes with id.
int n_values
Number of values in the list.
int * value
Array of values.
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.
void RTreeSetOverflow(struct RTree *t, char overflow)
Enable/disable R*-tree forced reinsertion (overflow)
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.
struct RTree * RTreeCreateTree(int fd, off_t rootpos, int ndims)
Create new empty R*-Tree.