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);
91 #define D ((ax2 - ax1) * (by1 - by2) - (ay2 - ay1) * (bx1 - bx2))
92 #define D1 ((bx1 - ax1) * (by1 - by2) - (by1 - ay1) * (bx1 - bx2))
93 #define D2 ((ax2 - ax1) * (by1 - ay1) - (ay2 - ay1) * (bx1 - ax1))
112 double ay2,
double az2,
double bx1,
double by1,
113 double bz1,
double bx2,
double by2,
double bz2,
114 double *x1,
double *y1,
double *z1,
double *x2,
115 double *y2,
double *z2,
int with_z)
117 static int first_3d = 1;
118 double d, d1, d2, r1, dtol,
t;
124 G_debug(4,
"Vect_segment_intersection()");
125 G_debug(4,
" %.15g , %.15g - %.15g , %.15g", ax1, ay1, ax2, ay2);
126 G_debug(4,
" %.15g , %.15g - %.15g , %.15g", bx1, by1, bx2, by2);
129 if (with_z && first_3d) {
130 G_warning(
_(
"3D not supported by Vect_segment_intersection()"));
146 else if (bx2 == bx1) {
166 else if (ax2 == ax1) {
184 if (ax1 == bx1 && ay1 == by1 && ax2 == bx2 && ay2 == by2) {
185 G_debug(2,
" -> identical segments");
200 else if (bx1 == ax1) {
203 else if (bx2 == ax2) {
206 else if (by1 == ay1) {
240 G_debug(2,
"Vect_segment_intersection(): d = %f, d1 = %f, d2 = %f", d, d1,
244 if (ax1 == bx1 && ay1 == by1) {
249 if (ax1 == bx2 && ay1 == by2) {
254 if (ax2 == bx1 && ay2 == by1) {
259 if (ax2 == bx2 && ay2 == by2) {
268 if (fabs(d) > dtol) {
270 G_debug(2,
" -> not parallel/collinear: d1 = %f, d2 = %f", d1, d2);
272 if (d1 < 0 || d1 > d || d2 < 0 || d2 > d) {
276 " -> fp error, but intersection at end points %f, %f",
282 G_debug(2,
" -> no intersection");
289 if (d1 < d || d1 > 0 || d2 < d || d2 > 0) {
293 " -> fp error, but intersection at end points %f, %f",
299 G_debug(2,
" -> no intersection");
308 *x1 = ax1 + r1 * (ax2 - ax1);
309 *y1 = ay1 + r1 * (ay2 - ay1);
312 G_debug(2,
" -> intersection %f, %f", *x1, *y1);
317 G_debug(3,
" -> parallel/collinear");
322 G_debug(2,
"Segments are apparently parallel, but connected at end "
323 "points -> collinear");
334 G_debug(2,
" -> collinear vertical");
335 if (ay1 > by2 || ay2 < by1) {
336 G_debug(2,
" -> no intersection");
342 G_debug(2,
" -> connected by end points");
350 G_debug(2,
" -> connected by end points");
359 G_debug(3,
" -> vertical overlap");
361 if (ay1 <= by1 && ay2 >= by2) {
362 G_debug(2,
" -> a contains b");
373 if (ay1 >= by1 && ay2 <= by2) {
374 G_debug(2,
" -> b contains a");
386 G_debug(2,
" -> partial overlap");
387 if (by1 > ay1 && by1 < ay2) {
398 if (by2 > ay1 && by2 < ay2) {
412 "Vect_segment_intersection() ERROR (collinear vertical segments)"));
425 G_debug(2,
" -> collinear non vertical");
428 if ((bx1 > ax2) || (bx2 < ax1)) {
430 G_debug(2,
" -> no intersection");
435 G_debug(2,
" -> overlap/connected end points");
438 if (ax1 == bx2 && ay1 == by2) {
439 G_debug(2,
" -> connected by end points");
446 if (ax2 == bx1 && ay2 == by1) {
447 G_debug(2,
" -> connected by end points");
456 if (ax1 <= bx1 && ax2 >= bx2) {
457 G_debug(2,
" -> a contains b");
468 if (ax1 >= bx1 && ax2 <= bx2) {
469 G_debug(2,
" -> b contains a");
481 G_debug(2,
" -> partial overlap");
482 if (bx1 > ax1 && bx1 < ax2) {
493 if (bx2 > ax1 && bx2 < ax2) {
507 "Vect_segment_intersection() ERROR (collinear non vertical segments)"));
528 static int a_cross = 0;
530 static CROSS *cross =
NULL;
531 static int *use_cross =
NULL;
533 static void add_cross(
int asegment,
double adistance,
int bsegment,
534 double bdistance,
double x,
double y)
536 if (n_cross == a_cross) {
539 (CROSS *)
G_realloc((
void *)cross, (a_cross + 101) *
sizeof(CROSS));
541 (
int *)
G_realloc((
void *)use_cross, (a_cross + 101) *
sizeof(
int));
547 " add new cross: aseg/dist = %d/%f bseg/dist = %d/%f, x = %f y = %f",
548 asegment, adistance, bsegment, bdistance,
x, y);
549 cross[n_cross].segment[0] = asegment;
550 cross[n_cross].distance[0] = adistance;
551 cross[n_cross].segment[1] = bsegment;
552 cross[n_cross].distance[1] = bdistance;
553 cross[n_cross].x =
x;
554 cross[n_cross].y = y;
558 static int cmp_cross(
const void *pa,
const void *pb)
560 CROSS *p1 = (CROSS *)pa;
561 CROSS *p2 = (CROSS *)pb;
563 if (p1->segment[current] < p2->segment[current])
565 if (p1->segment[current] > p2->segment[current])
568 if (p1->distance[current] < p2->distance[current])
570 if (p1->distance[current] > p2->distance[current])
575 static double dist2(
double x1,
double y1,
double x2,
double y2)
581 return (dx * dx + dy * dy);
586 static int ident(
double x1,
double y1,
double x2,
double y2,
double thresh)
592 if ((dx * dx + dy * dy) <= thresh * thresh)
604 static int cross_seg(
int id,
const struct RTree_Rect *rect
UNUSED,
void *arg)
606 double x1, y1, z1, x2, y2, z2;
615 APnts->
x[i], APnts->
y[i], APnts->
z[i], APnts->
x[i + 1], APnts->
y[i + 1],
616 APnts->
z[i + 1], BPnts->
x[j], BPnts->
y[j], BPnts->
z[j], BPnts->
x[j + 1],
617 BPnts->
y[j + 1], BPnts->
z[j + 1], &x1, &y1, &z1, &x2, &y2, &z2, 0);
621 G_debug(2,
" -> %d x %d: intersection type = %d", i, j, ret);
623 G_debug(3,
" in %f, %f ", x1, y1);
624 add_cross(i, 0.0, j, 0.0, x1, y1);
626 else if (ret == 2 || ret == 3 || ret == 4 || ret == 5) {
631 G_debug(3,
" in %f, %f; %f, %f", x1, y1, x2, y2);
632 add_cross(i, 0.0, j, 0.0, x1, y1);
633 add_cross(i, 0.0, j, 0.0, x2, y2);
665 struct line_pnts ***BLines,
int *nalines,
666 int *nblines,
int with_z
UNUSED)
668 int i, j, k,
l, last_seg, seg, last;
670 double dist, curdist, last_x, last_y, last_z;
671 double x,
y, rethresh;
673 struct RTree *MyRTree;
675 int seg1, seg2, vert1, vert2;
677 static int rect_init = 0;
680 if (debug_level == -1) {
684 debug_level = atoi(dstr);
768 if (abbox.
N > ABox->
N)
770 if (abbox.
S < ABox->
S)
772 if (abbox.
E > ABox->
E)
774 if (abbox.
W < ABox->
W)
776 if (abbox.
T > ABox->
T)
778 if (abbox.
B < ABox->
B)
791 for (i = 0; i < BPoints->
n_points - 1; i++) {
792 if (BPoints->
x[i] <= BPoints->
x[i + 1]) {
801 if (BPoints->
y[i] <= BPoints->
y[i + 1]) {
810 if (BPoints->
z[i] <= BPoints->
z[i + 1]) {
834 for (i = 0; i < APoints->
n_points - 1; i++) {
835 if (APoints->
x[i] <= APoints->
x[i + 1]) {
844 if (APoints->
y[i] <= APoints->
y[i + 1]) {
852 if (APoints->
z[i] <= APoints->
z[i + 1]) {
876 G_debug(2,
"n_cross = %d", n_cross);
886 for (i = 0; i < n_cross; i++) {
889 seg = cross[i].segment[0];
891 dist2(cross[i].
x, cross[i].y, APoints->
x[seg], APoints->
y[seg]);
895 cross[i].distance[0] = curdist;
898 dist = dist2(cross[i].
x, cross[i].y, APoints->
x[seg + 1],
899 APoints->
y[seg + 1]);
900 if (dist < curdist) {
902 x = APoints->
x[seg + 1];
903 y = APoints->
y[seg + 1];
907 seg = cross[i].segment[1];
908 dist = dist2(cross[i].
x, cross[i].y, BPoints->
x[seg], BPoints->
y[seg]);
909 cross[i].distance[1] = dist;
911 if (dist < curdist) {
917 dist = dist2(cross[i].
x, cross[i].y, BPoints->
x[seg + 1],
918 BPoints->
y[seg + 1]);
919 if (dist < curdist) {
921 x = BPoints->
x[seg + 1];
922 y = BPoints->
y[seg + 1];
924 if (curdist < rethresh * rethresh) {
929 seg = cross[i].segment[0];
930 cross[i].distance[0] =
931 dist2(APoints->
x[seg], APoints->
y[seg], cross[i].x, cross[i].y);
932 seg = cross[i].segment[1];
933 cross[i].distance[1] =
934 dist2(BPoints->
x[seg], BPoints->
y[seg], cross[i].x, cross[i].y);
939 for (
l = 1;
l < 3;
l++) {
940 for (i = 0; i < n_cross; i++)
947 G_debug(2,
"Clean and create array for line A");
955 G_debug(2,
"Clean and create array for line B");
964 qsort((
void *)cross,
sizeof(
char) * n_cross,
sizeof(CROSS), cmp_cross);
968 if (debug_level > 2) {
969 for (i = 0; i < n_cross; i++) {
972 " cross = %d seg1/dist1 = %d/%f seg2/dist2 = %d/%f x = %f "
974 i, cross[i].segment[current],
975 sqrt(cross[i].distance[current]), cross[i].segment[second],
976 sqrt(cross[i].distance[second]), cross[i].
x, cross[i].y);
981 for (i = 0; i < n_cross; i++) {
982 if (use_cross[i] == 1) {
986 if ((cross[i].segment[current] == 0 &&
987 cross[i].
x == Points1->
x[0] &&
988 cross[i].y == Points1->
y[0]) ||
989 (cross[i].segment[current] == j - 1 &&
990 cross[i].x == Points1->
x[j] &&
991 cross[i].y == Points1->
y[j])) {
993 G_debug(3,
"cross %d deleted (first/last point)", i);
1010 for (i = 0; i < n_cross; i++) {
1011 if (use_cross[i] == 0)
1013 G_debug(3,
" is %d between colinear?", i);
1015 seg1 = cross[i].segment[current];
1016 seg2 = cross[i].segment[second];
1019 if (cross[i].
x == Points1->
x[seg1] &&
1020 cross[i].y == Points1->
y[seg1]) {
1023 else if (cross[i].
x == Points1->
x[seg1 + 1] &&
1024 cross[i].y == Points1->
y[seg1 + 1]) {
1028 G_debug(3,
" -> is not vertex on 1. line");
1036 if (cross[i].
x == Points2->
x[seg2] &&
1037 cross[i].y == Points2->
y[seg2]) {
1040 else if (cross[i].
x == Points2->
x[seg2 + 1] &&
1041 cross[i].y == Points2->
y[seg2 + 1]) {
1045 G_debug(3,
" -> is not vertex on 2. line");
1048 G_debug(3,
" seg1/vert1 = %d/%d seg2/ver2 = %d/%d", seg1, vert1,
1052 if (vert2 == 0 || vert2 == Points2->
n_points - 1) {
1053 G_debug(3,
" -> vertex 2 (%d) is first/last", vert2);
1058 if (!((Points1->
x[vert1 - 1] == Points2->
x[vert2 - 1] &&
1059 Points1->
y[vert1 - 1] == Points2->
y[vert2 - 1] &&
1060 Points1->
x[vert1 + 1] == Points2->
x[vert2 + 1] &&
1061 Points1->
y[vert1 + 1] == Points2->
y[vert2 + 1]) ||
1062 (Points1->
x[vert1 - 1] == Points2->
x[vert2 + 1] &&
1063 Points1->
y[vert1 - 1] == Points2->
y[vert2 + 1] &&
1064 Points1->
x[vert1 + 1] == Points2->
x[vert2 - 1] &&
1065 Points1->
y[vert1 + 1] == Points2->
y[vert2 - 1]))) {
1066 G_debug(3,
" -> previous/next are not identical");
1072 G_debug(3,
" -> collinear -> remove");
1093 for (i = 0; i < n_cross; i++) {
1094 if (use_cross[i] == 0)
1100 seg = cross[i].segment[current];
1102 G_debug(3,
" duplicate ?: cross = %d seg = %d dist = %f", i,
1103 cross[i].segment[current], cross[i].distance[current]);
1104 if ((cross[i].segment[current] == cross[last].segment[current] &&
1105 cross[i].distance[current] == cross[last].distance[current]) ||
1106 (cross[i].segment[current] ==
1107 cross[last].segment[current] + 1 &&
1108 cross[i].distance[current] == 0 &&
1109 cross[i].
x == cross[last].
x && cross[i].y == cross[last].y)) {
1110 G_debug(3,
" cross %d identical to last -> removed", i);
1121 G_debug(3,
" alive crosses:");
1122 for (i = 0; i < n_cross; i++) {
1123 if (use_cross[i] == 1) {
1130 if (n_alive_cross > 0) {
1132 use_cross[n_cross] = 1;
1134 cross[n_cross].x = Points->
x[j];
1135 cross[n_cross].y = Points->
y[j];
1136 cross[n_cross].segment[current] = Points->
n_points - 2;
1139 last_x = Points->
x[0];
1140 last_y = Points->
y[0];
1141 last_z = Points->
z[0];
1144 for (i = 0; i <= n_cross; i++) {
1145 seg = cross[i].segment[current];
1146 G_debug(2,
"%d seg = %d dist = %f", i, seg,
1147 cross[i].distance[current]);
1148 if (use_cross[i] == 0) {
1149 G_debug(3,
" removed -> next");
1157 G_debug(2,
" append last vert: %f %f", last_x, last_y);
1160 for (j = last_seg + 1; j <= seg; j++) {
1161 G_debug(2,
" segment j = %d", j);
1163 if ((j == last_seg + 1) && Points->
x[j] == last_x &&
1164 Points->
y[j] == last_y) {
1165 G_debug(2,
" -> skip (identical to last break)");
1170 G_debug(2,
" append first of seg: %f %f", Points->
x[j],
1175 last_x = cross[i].x;
1176 last_y = cross[i].y;
1179 if (Points->
z[last_seg] == Points->
z[last_seg + 1]) {
1180 last_z = Points->
z[last_seg + 1];
1182 else if (last_x == Points->
x[last_seg] &&
1183 last_y == Points->
y[last_seg]) {
1184 last_z = Points->
z[last_seg];
1186 else if (last_x == Points->
x[last_seg + 1] &&
1187 last_y == Points->
y[last_seg + 1]) {
1188 last_z = Points->
z[last_seg + 1];
1191 dist = dist2(Points->
x[last_seg], Points->
x[last_seg + 1],
1192 Points->
y[last_seg], Points->
y[last_seg + 1]);
1194 (Points->
z[last_seg] *
1195 sqrt(cross[i].distance[current]) +
1196 Points->
z[last_seg + 1] *
1197 (sqrt(dist) - sqrt(cross[i].distance[current]))) /
1203 G_debug(2,
" append cross / last point: %f %f", cross[i].
x,
1208 G_debug(2,
" line is degenerate -> skipped");
1231 static struct line_pnts *APnts, *BPnts, *IPnts;
1233 static int cross_found;
1234 static int report_all;
1237 static int find_cross(
int id,
const struct RTree_Rect *rect
UNUSED,
void *arg)
1239 double x1, y1, z1, x2, y2, z2;
1248 APnts->
x[i], APnts->
y[i], APnts->
z[i], APnts->
x[i + 1], APnts->
y[i + 1],
1249 APnts->
z[i + 1], BPnts->
x[j], BPnts->
y[j], BPnts->
z[j], BPnts->
x[j + 1],
1250 BPnts->
y[j + 1], BPnts->
z[j + 1], &x1, &y1, &z1, &x2, &y2, &z2, 0);
1258 G_warning(
_(
"Error while adding point to array. Out of memory"));
1264 G_warning(
_(
"Error while adding point to array. Out of memory"));
1266 G_warning(
_(
"Error while adding point to array. Out of memory"));
1282 double dist, rethresh;
1283 struct RTree *MyRTree;
1285 static int rect_init = 0;
1292 rethresh = 0.000001;
1305 if (APoints->
x[0] == BPoints->
x[0] && APoints->
y[0] == BPoints->
y[0]) {
1309 &APoints->
y[0],
NULL, 1))
1311 _(
"Error while adding point to array. Out of memory"));
1315 if (APoints->
z[0] == BPoints->
z[0]) {
1320 G_warning(
_(
"Error while adding point to array. Out of "
1337 if (dist <= rethresh) {
1342 _(
"Error while adding point to array. Out of memory"));
1354 if (dist <= rethresh) {
1359 _(
"Error while adding point to array. Out of memory"));
1376 for (i = 0; i < BPoints->
n_points - 1; i++) {
1377 if (BPoints->
x[i] <= BPoints->
x[i + 1]) {
1386 if (BPoints->
y[i] <= BPoints->
y[i + 1]) {
1395 if (BPoints->
z[i] <= BPoints->
z[i + 1]) {
1412 for (i = 0; i < APoints->
n_points - 1; i++) {
1413 if (APoints->
x[i] <= APoints->
x[i + 1]) {
1422 if (APoints->
y[i] <= APoints->
y[i + 1]) {
1430 if (APoints->
z[i] <= APoints->
z[i + 1]) {
1442 if (!report_all && cross_found) {
void G_warning(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
const char * G_getenv_nofatal(const char *)
Get environment variable.
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
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_distance(const struct line_pnts *, double, double, double, int, double *, double *, double *, double *, double *, double *)
Calculate distance of point to line.
int Vect_box_overlap(const struct bound_box *, const struct bound_box *)
Tests for overlap of two boxes.
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_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
int dig_line_degenerate(const struct line_pnts *)
#define UNUSED
A macro for an attribute, if attached to a variable, indicating that the variable is not used.
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.
int line_check_intersection(struct line_pnts *APoints, struct line_pnts *BPoints, int with_z)
int Vect_line_check_intersection(struct line_pnts *APoints, struct line_pnts *BPoints, int with_z)
Check if 2 lines intersect.
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.
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 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 UNUSED)
Intersect 2 lines.
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.