76 static int cmp_cross(
const void *pa,
const void *pb);
77 static void add_cross(
int asegment,
double adistance,
int bsegment,
78 double bdistance,
double x,
double y);
79 static double dist2(
double x1,
double y1,
double x2,
double y2);
81 static int debug_level = -1;
84 static int ident(
double x1,
double y1,
double x2,
double y2,
double thresh);
86 static int cross_seg(
int id,
const struct RTree_Rect *rect,
void *arg);
87 static int find_cross(
int id,
const struct RTree_Rect *rect,
void *arg);
89 #define D ((ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2)) 90 #define D1 ((bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2)) 91 #define D2 ((ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1)) 110 double ay2,
double az2,
double bx1,
double by1,
111 double bz1,
double bx2,
double by2,
double bz2,
112 double *x1,
double *y1,
double *z1,
double *x2,
113 double *y2,
double *z2,
int with_z)
115 static int first_3d = 1;
116 double d, d1, d2, r1, dtol,
t;
122 G_debug(4,
"Vect_segment_intersection()");
123 G_debug(4,
" %.15g , %.15g - %.15g , %.15g", ax1, ay1, ax2, ay2);
124 G_debug(4,
" %.15g , %.15g - %.15g , %.15g", bx1, by1, bx2, by2);
127 if (with_z && first_3d) {
128 G_warning(
_(
"3D not supported by Vect_segment_intersection()"));
144 else if (bx2 == bx1) {
164 else if (ax2 == ax1) {
182 if (ax1 == bx1 && ay1 == by1 && ax2 == bx2 && ay2 == by2) {
183 G_debug(2,
" -> identical segments");
198 else if (bx1 == ax1) {
201 else if (bx2 == ax2) {
204 else if (by1 == ay1) {
238 G_debug(2,
"Vect_segment_intersection(): d = %f, d1 = %f, d2 = %f", d, d1,
242 if (ax1 == bx1 && ay1 == by1) {
247 if (ax1 == bx2 && ay1 == by2) {
252 if (ax2 == bx1 && ay2 == by1) {
257 if (ax2 == bx2 && ay2 == by2) {
266 if (fabs(d) > dtol) {
268 G_debug(2,
" -> not parallel/collinear: d1 = %f, d2 = %f", d1, d2);
270 if (d1 < 0 || d1 > d || d2 < 0 || d2 > d) {
272 G_debug(2,
" -> fp error, but intersection at end points %f, %f", *x1, *y1);
277 G_debug(2,
" -> no intersection");
284 if (d1 < d || d1 > 0 || d2 < d || d2 > 0) {
286 G_debug(2,
" -> fp error, but intersection at end points %f, %f", *x1, *y1);
291 G_debug(2,
" -> no intersection");
300 *x1 = ax1 + r1 * (ax2 - ax1);
301 *y1 = ay1 + r1 * (ay2 - ay1);
304 G_debug(2,
" -> intersection %f, %f", *x1, *y1);
309 G_debug(3,
" -> parallel/collinear");
314 G_debug(2,
"Segments are apparently parallel, but connected at end points -> collinear");
325 G_debug(2,
" -> collinear vertical");
326 if (ay1 > by2 || ay2 < by1) {
327 G_debug(2,
" -> no intersection");
333 G_debug(2,
" -> connected by end points");
341 G_debug(2,
" -> connected by end points");
350 G_debug(3,
" -> vertical overlap");
352 if (ay1 <= by1 && ay2 >= by2) {
353 G_debug(2,
" -> a contains b");
364 if (ay1 >= by1 && ay2 <= by2) {
365 G_debug(2,
" -> b contains a");
377 G_debug(2,
" -> partial overlap");
378 if (by1 > ay1 && by1 < ay2) {
389 if (by2 > ay1 && by2 < ay2) {
402 G_warning(
_(
"Vect_segment_intersection() ERROR (collinear vertical segments)"));
415 G_debug(2,
" -> collinear non vertical");
418 if ((bx1 > ax2) || (bx2 < ax1)) {
420 G_debug(2,
" -> no intersection");
425 G_debug(2,
" -> overlap/connected end points");
428 if (ax1 == bx2 && ay1 == by2) {
429 G_debug(2,
" -> connected by end points");
436 if (ax2 == bx1 && ay2 == by1) {
437 G_debug(2,
" -> connected by end points");
446 if (ax1 <= bx1 && ax2 >= bx2) {
447 G_debug(2,
" -> a contains b");
458 if (ax1 >= bx1 && ax2 <= bx2) {
459 G_debug(2,
" -> b contains a");
471 G_debug(2,
" -> partial overlap");
472 if (bx1 > ax1 && bx1 < ax2) {
483 if (bx2 > ax1 && bx2 < ax2) {
496 G_warning(
_(
"Vect_segment_intersection() ERROR (collinear non vertical segments)"));
518 static int a_cross = 0;
520 static CROSS *cross =
NULL;
521 static int *use_cross =
NULL;
523 static void add_cross(
int asegment,
double adistance,
int bsegment,
524 double bdistance,
double x,
double y)
526 if (n_cross == a_cross) {
530 (a_cross + 101) *
sizeof(CROSS));
533 (a_cross + 101) *
sizeof(
int));
538 " add new cross: aseg/dist = %d/%f bseg/dist = %d/%f, x = %f y = %f",
539 asegment, adistance, bsegment, bdistance,
x, y);
540 cross[n_cross].segment[0] = asegment;
541 cross[n_cross].distance[0] = adistance;
542 cross[n_cross].segment[1] = bsegment;
543 cross[n_cross].distance[1] = bdistance;
544 cross[n_cross].x =
x;
545 cross[n_cross].y = y;
549 static int cmp_cross(
const void *pa,
const void *pb)
551 CROSS *p1 = (CROSS *) pa;
552 CROSS *p2 = (CROSS *) pb;
554 if (p1->segment[current] < p2->segment[current])
556 if (p1->segment[current] > p2->segment[current])
559 if (p1->distance[current] < p2->distance[current])
561 if (p1->distance[current] > p2->distance[current])
566 static double dist2(
double x1,
double y1,
double x2,
double y2)
572 return (dx * dx + dy * dy);
577 static int ident(
double x1,
double y1,
double x2,
double y2,
double thresh)
583 if ((dx * dx + dy * dy) <= thresh * thresh)
594 static int cross_seg(
int id,
const struct RTree_Rect *rect,
void *arg)
596 double x1, y1, z1, x2, y2, z2;
605 APnts->
x[i + 1], APnts->
y[i + 1],
606 APnts->
z[i + 1], BPnts->
x[j], BPnts->
y[j],
607 BPnts->
z[j], BPnts->
x[j + 1],
608 BPnts->
y[j + 1], BPnts->
z[j + 1], &x1,
609 &y1, &z1, &x2, &y2, &z2, 0);
613 G_debug(2,
" -> %d x %d: intersection type = %d", i, j, ret);
615 G_debug(3,
" in %f, %f ", x1, y1);
616 add_cross(i, 0.0, j, 0.0, x1, y1);
618 else if (ret == 2 || ret == 3 || ret == 4 || ret == 5) {
623 G_debug(3,
" in %f, %f; %f, %f", x1, y1, x2, y2);
624 add_cross(i, 0.0, j, 0.0, x1, y1);
625 add_cross(i, 0.0, j, 0.0, x2, y2);
656 int *nalines,
int *nblines,
int with_z)
658 int i, j, k,
l, last_seg, seg, last;
660 double dist, curdist, last_x, last_y, last_z;
661 double x,
y, rethresh;
663 struct RTree *MyRTree;
665 int seg1, seg2, vert1, vert2;
667 static int rect_init = 0;
670 if (debug_level == -1) {
674 debug_level = atoi(dstr);
756 if (abbox.
N > ABox->
N)
758 if (abbox.
S < ABox->
S)
760 if (abbox.
E > ABox->
E)
762 if (abbox.
W < ABox->
W)
764 if (abbox.
T > ABox->
T)
766 if (abbox.
B < ABox->
B)
779 for (i = 0; i < BPoints->
n_points - 1; i++) {
780 if (BPoints->
x[i] <= BPoints->
x[i + 1]) {
789 if (BPoints->
y[i] <= BPoints->
y[i + 1]) {
798 if (BPoints->
z[i] <= BPoints->
z[i + 1]) {
820 for (i = 0; i < APoints->
n_points - 1; i++) {
821 if (APoints->
x[i] <= APoints->
x[i + 1]) {
830 if (APoints->
y[i] <= APoints->
y[i + 1]) {
838 if (APoints->
z[i] <= APoints->
z[i + 1]) {
861 G_debug(2,
"n_cross = %d", n_cross);
871 for (i = 0; i < n_cross; i++) {
874 seg = cross[i].segment[0];
876 dist2(cross[i].x, cross[i].y, APoints->
x[seg], APoints->
y[seg]);
880 cross[i].distance[0] = curdist;
884 dist2(cross[i].x, cross[i].y, APoints->
x[seg + 1],
885 APoints->
y[seg + 1]);
886 if (dist < curdist) {
888 x = APoints->
x[seg + 1];
889 y = APoints->
y[seg + 1];
893 seg = cross[i].segment[1];
895 dist2(cross[i].x, cross[i].y, BPoints->
x[seg], BPoints->
y[seg]);
896 cross[i].distance[1] = dist;
898 if (dist < curdist) {
904 dist = dist2(cross[i].x, cross[i].y, BPoints->
x[seg + 1], BPoints->
y[seg + 1]);
905 if (dist < curdist) {
907 x = BPoints->
x[seg + 1];
908 y = BPoints->
y[seg + 1];
910 if (curdist < rethresh * rethresh) {
915 seg = cross[i].segment[0];
916 cross[i].distance[0] =
917 dist2(APoints->
x[seg], APoints->
y[seg], cross[i].x, cross[i].y);
918 seg = cross[i].segment[1];
919 cross[i].distance[1] =
920 dist2(BPoints->
x[seg], BPoints->
y[seg], cross[i].x, cross[i].y);
925 for (l = 1; l < 3; l++) {
926 for (i = 0; i < n_cross; i++)
933 G_debug(2,
"Clean and create array for line A");
941 G_debug(2,
"Clean and create array for line B");
950 qsort((
void *)cross,
sizeof(
char) * n_cross,
sizeof(CROSS),
955 if (debug_level > 2) {
956 for (i = 0; i < n_cross; i++) {
958 " cross = %d seg1/dist1 = %d/%f seg2/dist2 = %d/%f x = %f y = %f",
959 i, cross[i].segment[current],
960 sqrt(cross[i].distance[current]),
961 cross[i].segment[second], sqrt(cross[i].distance[second]),
962 cross[i].x, cross[i].y);
967 for (i = 0; i < n_cross; i++) {
968 if (use_cross[i] == 1) {
972 if ((cross[i].segment[current] == 0 &&
973 cross[i].x == Points1->
x[0] &&
974 cross[i].y == Points1->
y[0]) ||
975 (cross[i].segment[current] == j - 1 &&
976 cross[i].x == Points1->
x[j] &&
977 cross[i].y == Points1->
y[j])) {
979 G_debug(3,
"cross %d deleted (first/last point)", i);
995 for (i = 0; i < n_cross; i++) {
996 if (use_cross[i] == 0)
998 G_debug(3,
" is %d between colinear?", i);
1000 seg1 = cross[i].segment[current];
1001 seg2 = cross[i].segment[second];
1004 if (cross[i].x == Points1->
x[seg1] &&
1005 cross[i].y == Points1->
y[seg1]) {
1008 else if (cross[i].x == Points1->
x[seg1 + 1] &&
1009 cross[i].y == Points1->
y[seg1 + 1]) {
1013 G_debug(3,
" -> is not vertex on 1. line");
1020 if (cross[i].x == Points2->
x[seg2] &&
1021 cross[i].y == Points2->
y[seg2]) {
1024 else if (cross[i].x == Points2->
x[seg2 + 1] &&
1025 cross[i].y == Points2->
y[seg2 + 1]) {
1029 G_debug(3,
" -> is not vertex on 2. line");
1032 G_debug(3,
" seg1/vert1 = %d/%d seg2/ver2 = %d/%d", seg1,
1033 vert1, seg2, vert2);
1036 if (vert2 == 0 || vert2 == Points2->
n_points - 1) {
1037 G_debug(3,
" -> vertex 2 (%d) is first/last", vert2);
1042 if (!((Points1->
x[vert1 - 1] == Points2->
x[vert2 - 1] &&
1043 Points1->
y[vert1 - 1] == Points2->
y[vert2 - 1] &&
1044 Points1->
x[vert1 + 1] == Points2->
x[vert2 + 1] &&
1045 Points1->
y[vert1 + 1] == Points2->
y[vert2 + 1]) ||
1046 (Points1->
x[vert1 - 1] == Points2->
x[vert2 + 1] &&
1047 Points1->
y[vert1 - 1] == Points2->
y[vert2 + 1] &&
1048 Points1->
x[vert1 + 1] == Points2->
x[vert2 - 1] &&
1049 Points1->
y[vert1 + 1] == Points2->
y[vert2 - 1])
1052 G_debug(3,
" -> previous/next are not identical");
1058 G_debug(3,
" -> collinear -> remove");
1077 for (i = 0; i < n_cross; i++) {
1078 if (use_cross[i] == 0)
1084 seg = cross[i].segment[current];
1086 G_debug(3,
" duplicate ?: cross = %d seg = %d dist = %f", i,
1087 cross[i].segment[current], cross[i].distance[current]);
1088 if ((cross[i].segment[current] == cross[last].segment[current] &&
1089 cross[i].distance[current] == cross[last].distance[current])
1090 || (cross[i].segment[current] ==
1091 cross[last].segment[current] + 1 &&
1092 cross[i].distance[current] == 0 &&
1093 cross[i].x == cross[last].x &&
1094 cross[i].y == cross[last].y)) {
1095 G_debug(3,
" cross %d identical to last -> removed", i);
1106 G_debug(3,
" alive crosses:");
1107 for (i = 0; i < n_cross; i++) {
1108 if (use_cross[i] == 1) {
1115 if (n_alive_cross > 0) {
1117 use_cross[n_cross] = 1;
1119 cross[n_cross].x = Points->
x[j];
1120 cross[n_cross].y = Points->
y[j];
1121 cross[n_cross].segment[current] = Points->
n_points - 2;
1124 last_x = Points->
x[0];
1125 last_y = Points->
y[0];
1126 last_z = Points->
z[0];
1129 for (i = 0; i <= n_cross; i++) {
1130 seg = cross[i].segment[current];
1131 G_debug(2,
"%d seg = %d dist = %f", i, seg,
1132 cross[i].distance[current]);
1133 if (use_cross[i] == 0) {
1134 G_debug(3,
" removed -> next");
1142 G_debug(2,
" append last vert: %f %f", last_x, last_y);
1145 for (j = last_seg + 1; j <= seg; j++) {
1146 G_debug(2,
" segment j = %d", j);
1148 if ((j == last_seg + 1) && Points->
x[j] == last_x &&
1149 Points->
y[j] == last_y) {
1150 G_debug(2,
" -> skip (identical to last break)");
1155 G_debug(2,
" append first of seg: %f %f", Points->
x[j],
1160 last_x = cross[i].x;
1161 last_y = cross[i].y;
1164 if (Points->
z[last_seg] == Points->
z[last_seg + 1]) {
1165 last_z = Points->
z[last_seg + 1];
1167 else if (last_x == Points->
x[last_seg] &&
1168 last_y == Points->
y[last_seg]) {
1169 last_z = Points->
z[last_seg];
1171 else if (last_x == Points->
x[last_seg + 1] &&
1172 last_y == Points->
y[last_seg + 1]) {
1173 last_z = Points->
z[last_seg + 1];
1176 dist = dist2(Points->
x[last_seg],
1177 Points->
x[last_seg + 1],
1178 Points->
y[last_seg],
1179 Points->
y[last_seg + 1]);
1180 last_z = (Points->
z[last_seg] * sqrt(cross[i].distance[current]) +
1181 Points->
z[last_seg + 1] * (sqrt(dist) - sqrt(cross[i].distance[current]))) /
1187 G_debug(2,
" append cross / last point: %f %f", cross[i].x,
1192 G_debug(2,
" line is degenerate -> skipped");
1215 static struct line_pnts *APnts, *BPnts, *IPnts;
1217 static int cross_found;
1220 static int find_cross(
int id,
const struct RTree_Rect *rect,
void *arg)
1222 double x1, y1, z1, x2, y2, z2;
1231 APnts->
x[i + 1], APnts->
y[i + 1],
1232 APnts->
z[i + 1], BPnts->
x[j], BPnts->
y[j],
1233 BPnts->
z[j], BPnts->
x[j + 1],
1234 BPnts->
y[j + 1], BPnts->
z[j + 1], &x1,
1235 &y1, &z1, &x2, &y2, &z2, 0);
1243 G_warning(
_(
"Error while adding point to array. Out of memory"));
1249 G_warning(
_(
"Error while adding point to array. Out of memory"));
1251 G_warning(
_(
"Error while adding point to array. Out of memory"));
1279 double dist, rethresh;
1280 struct RTree *MyRTree;
1282 static int rect_init = 0;
1289 rethresh = 0.000001;
1301 if (APoints->
x[0] == BPoints->
x[0] && APoints->
y[0] == BPoints->
y[0]) {
1305 &APoints->
y[0],
NULL, 1))
1306 G_warning(
_(
"Error while adding point to array. Out of memory"));
1310 if (APoints->
z[0] == BPoints->
z[0]) {
1313 &APoints->
y[0], &APoints->
z[0],
1315 G_warning(
_(
"Error while adding point to array. Out of memory"));
1332 if (dist <= rethresh) {
1336 G_warning(
_(
"Error while adding point to array. Out of memory"));
1349 if (dist <= rethresh) {
1353 G_warning(
_(
"Error while adding point to array. Out of memory"));
1370 for (i = 0; i < BPoints->
n_points - 1; i++) {
1371 if (BPoints->
x[i] <= BPoints->
x[i + 1]) {
1380 if (BPoints->
y[i] <= BPoints->
y[i + 1]) {
1389 if (BPoints->
z[i] <= BPoints->
z[i + 1]) {
1404 for (i = 0; i < APoints->
n_points - 1; i++) {
1405 if (APoints->
x[i] <= APoints->
x[i + 1]) {
1414 if (APoints->
y[i] <= APoints->
y[i + 1]) {
1422 if (APoints->
z[i] <= APoints->
z[i + 1]) {
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 Vect_line_get_intersections(struct line_pnts *APoints, struct line_pnts *BPoints, struct line_pnts *IPoints, int with_z)
Get 2 lines intersection points.
int RTreeSearch(struct RTree *t, struct RTree_Rect *r, SearchHitCallback *shcb, void *cbarg)
Search an R*-Tree.
int n_points
Number of points.
int Vect_copy_xyz_to_pnts(struct line_pnts *, const double *, const double *, const double *, int)
Copy points from array to line_pnts structure.
int Vect_line_check_intersection(struct line_pnts *APoints, struct line_pnts *BPoints, int with_z)
Check if 2 lines intersect.
int Vect_line_intersection(struct line_pnts *APoints, struct line_pnts *BPoints, struct bound_box *ABox, struct bound_box *BBox, struct line_pnts ***ALines, struct line_pnts ***BLines, int *nalines, int *nblines, int with_z)
Intersect 2 lines.
int Vect_line_distance(const struct line_pnts *, double, double, double, int, double *, double *, double *, double *, double *, double *)
Calculate distance of point to line.
double * x
Array of X coordinates.
Feature geometry info - coordinates.
int Vect_segment_intersection(double ax1, double ay1, double az1, double ax2, double ay2, double az2, double bx1, double by1, double bz1, double bx2, double by2, double bz2, double *x1, double *y1, double *z1, double *x2, double *y2, double *z2, int with_z)
Check for intersect of 2 line segments.
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
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.
double * y
Array of Y coordinates.
void G_warning(const char *,...) __attribute__((format(printf
double * z
Array of Z coordinates.
const char * G_getenv_nofatal(const char *)
Get environment variable.
int dig_line_degenerate(const struct line_pnts *)
int Vect_box_overlap(const struct bound_box *, const struct bound_box *)
Tests for overlap of two boxes.
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
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.