78 static int cmp_cross(
const void *pa,
const void *pb);
79 static void add_cross(
int asegment,
double adistance,
int bsegment,
80 double bdistance,
double x,
double y);
81 static double dist2(
double x1,
double y1,
double x2,
double y2);
83 static int debug_level = -1;
86 static int ident(
double x1,
double y1,
double x2,
double y2,
double thresh);
88 static int snap_cross(
int asegment,
double *adistance,
int bsegment,
89 double *bdistance,
double *xc,
double *yc);
90 static int cross_seg(
int i,
int j,
int b);
91 static int find_cross(
int i,
int j,
int b);
105 static int a_cross = 0;
107 static CROSS *cross =
NULL;
108 static int *use_cross =
NULL;
110 static double rethresh = 0.000001;
112 static double d_ulp(
double a,
double b)
127 result = frexp(dmax, &exp);
129 result = ldexp(result, exp);
134 static void add_cross(
int asegment,
double adistance,
int bsegment,
135 double bdistance,
double x,
double y)
137 if (n_cross == a_cross) {
141 (a_cross + 101) *
sizeof(CROSS));
144 (a_cross + 101) *
sizeof(
int));
149 " add new cross: aseg/dist = %d/%f bseg/dist = %d/%f, x = %f y = %f",
150 asegment, adistance, bsegment, bdistance, x, y);
151 cross[n_cross].segment[0] = asegment;
152 cross[n_cross].distance[0] = adistance;
153 cross[n_cross].segment[1] = bsegment;
154 cross[n_cross].distance[1] = bdistance;
155 cross[n_cross].x =
x;
156 cross[n_cross].y = y;
160 static int cmp_cross(
const void *pa,
const void *pb)
162 CROSS *p1 = (CROSS *) pa;
163 CROSS *p2 = (CROSS *) pb;
165 if (p1->segment[current] < p2->segment[current])
167 if (p1->segment[current] > p2->segment[current])
170 if (p1->distance[current] < p2->distance[current])
172 if (p1->distance[current] > p2->distance[current])
177 static double dist2(
double x1,
double y1,
double x2,
double y2)
183 return (dx * dx + dy * dy);
188 static int ident(
double x1,
double y1,
double x2,
double y2,
double thresh)
194 if ((dx * dx + dy * dy) <= thresh * thresh)
202 static struct line_pnts *APnts, *BPnts, *ABPnts[2], *IPnts;
206 static int snap_cross(
int asegment,
double *adistance,
int bsegment,
207 double *bdistance,
double *xc,
double *yc)
211 double dist, curdist, dthresh;
215 curdist = dist2(*xc, *yc, APnts->
x[seg], APnts->
y[seg]);
219 *adistance = curdist;
222 dist = dist2(*xc, *yc, APnts->
x[seg + 1], APnts->
y[seg + 1]);
223 if (dist < curdist) {
225 x = APnts->
x[seg + 1];
226 y = APnts->
y[seg + 1];
231 dist = dist2(*xc, *yc, BPnts->
x[seg], BPnts->
y[seg]);
234 if (dist < curdist) {
240 dist = dist2(*xc, *yc, BPnts->
x[seg + 1], BPnts->
y[seg + 1]);
241 if (dist < curdist) {
243 x = BPnts->
x[seg + 1];
244 y = BPnts->
y[seg + 1];
254 dthresh = d_ulp(x, y);
255 if (curdist < dthresh * dthresh) {
261 *adistance = dist2(*xc, *yc, APnts->
x[seg], APnts->
y[seg]);
263 *bdistance = dist2(*xc, *yc, BPnts->
x[seg], BPnts->
y[seg]);
272 static int cross_seg(
int i,
int j,
int b)
274 double x1, y1, z1, x2, y2, z2;
275 double y1min, y1max, y2min, y2max;
280 y1max = APnts->
y[i + 1];
281 if (APnts->
y[i] > APnts->
y[i + 1]) {
282 y1min = APnts->
y[i + 1];
287 y2max = BPnts->
y[j + 1];
288 if (BPnts->
y[j] > BPnts->
y[j + 1]) {
289 y2min = BPnts->
y[j + 1];
293 if (y1min > y2max || y1max < y2min)
299 APnts->
x[i + 1], APnts->
y[i + 1],
301 BPnts->
x[j], BPnts->
y[j],
303 BPnts->
x[j + 1], BPnts->
y[j + 1],
305 &x1, &y1, &z1, &x2, &y2, &z2, 0);
310 BPnts->
x[j + 1], BPnts->
y[j + 1],
312 APnts->
x[i], APnts->
y[i],
314 APnts->
x[i + 1], APnts->
y[i + 1],
316 &x1, &y1, &z1, &x2, &y2, &z2, 0);
321 G_debug(2,
" -> %d x %d: intersection type = %d", i, j, ret);
323 G_debug(3,
" in %f, %f ", x1, y1);
325 snap_cross(i, &adist, j, &bdist, &x1, &y1);
326 add_cross(i, adist, j, bdist, x1, y1);
328 add_cross(j, bdist, i, adist, x1, y1);
330 else if (ret == 2 || ret == 3 || ret == 4 || ret == 5) {
335 G_debug(3,
" in %f, %f; %f, %f", x1, y1, x2, y2);
336 snap_cross(i, &adist, j, &bdist, &x1, &y1);
337 add_cross(i, adist, j, bdist, x1, y1);
339 add_cross(j, bdist, i, adist, x1, y1);
340 snap_cross(i, &adist, j, &bdist, &x2, &y2);
341 add_cross(i, adist, j, bdist, x2, y2);
343 add_cross(j, bdist, i, adist, x2, y2);
354 #define GET_PARENT(p, c) ((p) = (int) (((c) - 2) / 3 + 1)) 355 #define GET_CHILD(c, p) ((c) = (int) (((p) * 3) - 1)) 374 static int cmp_q_x(
struct qitem *a,
struct qitem *b)
376 double x1, y1, z1, x2, y2, z2;
378 x1 = ABPnts[a->l]->
x[a->p];
379 y1 = ABPnts[a->l]->
y[a->p];
380 z1 = ABPnts[a->l]->
z[a->p];
382 x2 = ABPnts[b->l]->
x[b->p];
383 y2 = ABPnts[b->l]->
y[b->p];
384 z2 = ABPnts[b->l]->
z[b->p];
414 static int sift_up(
struct boq *q,
int start)
416 register int parent, child;
427 if (cmp_q_x(&a, b)) {
429 q->i[child] = q->i[parent];
445 static int boq_add(
struct boq *q,
struct qitem *i)
447 if (q->count + 2 >= q->alloc) {
448 q->alloc = q->count + 100;
449 q->i =
G_realloc(q->i, q->alloc *
sizeof(
struct qitem));
451 q->i[q->count + 1] = *i;
452 sift_up(q, q->count + 1);
460 static int boq_drop(
struct boq *q,
struct qitem *qi)
462 register int child, childr, parent;
481 while (
GET_CHILD(child, parent) <= q->count) {
484 for (childr = child + 1; childr <= q->count && childr < i; childr++) {
493 q->i[parent] = q->i[child];
498 if (parent < q->
count) {
499 q->i[parent] = q->i[q->count];
513 static int cmp_t_y(
const void *aa,
const void *bb)
515 double x1, y1, z1, x2, y2, z2;
516 struct qitem *a = (
struct qitem *) aa;
517 struct qitem *b = (
struct qitem *) bb;
519 x1 = ABPnts[a->l]->
x[a->p];
520 y1 = ABPnts[a->l]->
y[a->p];
521 z1 = ABPnts[a->l]->
z[a->p];
523 x2 = ABPnts[b->l]->
x[b->p];
524 y2 = ABPnts[b->l]->
y[b->p];
525 z2 = ABPnts[b->l]->
z[b->p];
555 double x1, y1, z1, x2, y2, z2;
563 for (i = 0; i < Pnts->
n_points - 1; i++) {
572 if (x1 == x2 && y1 == y2 && (!with_z || z1 == z2))
599 box.
W -= d_ulp(box.
W, box.
W);
600 box.
S -= d_ulp(box.
S, box.
S);
601 box.
B -= d_ulp(box.
B, box.
B);
602 box.
E += d_ulp(box.
E, box.
E);
603 box.
N += d_ulp(box.
N, box.
N);
604 box.
T += d_ulp(box.
T, box.
T);
689 int *nalines,
int *nblines,
int with_z)
691 int i, j, k,
l, nl, last_seg, seg, last;
693 double dist, last_x, last_y, last_z;
696 int seg1, seg2, vert1, vert2;
699 struct qitem qi, *found;
704 if (debug_level == -1) {
708 debug_level = atoi(dstr);
799 if (abbox.
N > ABox.
N)
801 if (abbox.
S < ABox.
S)
803 if (abbox.
E > ABox.
E)
805 if (abbox.
W < ABox.
W)
809 if (abbox.
T > BBox.
T)
811 if (abbox.
B < BBox.
B)
816 abbox.
N += d_ulp(abbox.
N, abbox.
N);
817 abbox.
S -= d_ulp(abbox.
S, abbox.
S);
818 abbox.
E += d_ulp(abbox.
E, abbox.
E);
819 abbox.
W -= d_ulp(abbox.
W, abbox.
W);
821 abbox.
T += d_ulp(abbox.
T, abbox.
T);
822 abbox.
B -= d_ulp(abbox.
B, abbox.
B);
826 G_fatal_error(
"Intersection with points is not yet supported");
833 bo_queue.i =
G_malloc(bo_queue.alloc *
sizeof(
struct qitem));
836 boq_load(&bo_queue, APnts, &abbox, 0, with_z);
840 boq_load(&bo_queue, BPnts, &abbox, 1, with_z);
851 while (boq_drop(&bo_queue, &qi)) {
859 cross_seg(qi.s, found->s, 0);
869 cross_seg(found->s, qi.s, 1);
900 G_debug(2,
"n_cross = %d", n_cross);
910 for (l = 1; l < nl; l++) {
911 for (i = 0; i < n_cross; i++)
918 G_debug(2,
"Clean and create array for line A");
926 G_debug(2,
"Clean and create array for line B");
935 qsort((
void *)cross,
sizeof(
char) * n_cross,
sizeof(CROSS),
940 if (debug_level > 2) {
941 for (i = 0; i < n_cross; i++) {
943 " cross = %d seg1/dist1 = %d/%f seg2/dist2 = %d/%f x = %f y = %f",
944 i, cross[i].segment[current],
945 sqrt(cross[i].distance[current]),
946 cross[i].segment[second], sqrt(cross[i].distance[second]),
947 cross[i].x, cross[i].y);
952 for (i = 0; i < n_cross; i++) {
953 if (use_cross[i] == 1) {
957 if ((cross[i].segment[current] == 0 &&
958 cross[i].x == Points1->
x[0] &&
959 cross[i].y == Points1->
y[0]) ||
960 (cross[i].segment[current] == j - 1 &&
961 cross[i].x == Points1->
x[j] &&
962 cross[i].y == Points1->
y[j])) {
964 G_debug(3,
"cross %d deleted (first/last point)", i);
980 for (i = 0; i < n_cross; i++) {
981 if (use_cross[i] == 0)
983 G_debug(3,
" is %d between colinear?", i);
985 seg1 = cross[i].segment[current];
986 seg2 = cross[i].segment[second];
989 if (cross[i].x == Points1->
x[seg1] &&
990 cross[i].y == Points1->
y[seg1]) {
993 else if (cross[i].x == Points1->
x[seg1 + 1] &&
994 cross[i].y == Points1->
y[seg1 + 1]) {
998 G_debug(3,
" -> is not vertex on 1. line");
1005 if (cross[i].x == Points2->
x[seg2] &&
1006 cross[i].y == Points2->
y[seg2]) {
1009 else if (cross[i].x == Points2->
x[seg2 + 1] &&
1010 cross[i].y == Points2->
y[seg2 + 1]) {
1014 G_debug(3,
" -> is not vertex on 2. line");
1017 G_debug(3,
" seg1/vert1 = %d/%d seg2/ver2 = %d/%d", seg1,
1018 vert1, seg2, vert2);
1021 if (vert2 == 0 || vert2 == Points2->
n_points - 1) {
1022 G_debug(3,
" -> vertex 2 (%d) is first/last", vert2);
1027 if (!((Points1->
x[vert1 - 1] == Points2->
x[vert2 - 1] &&
1028 Points1->
y[vert1 - 1] == Points2->
y[vert2 - 1] &&
1029 Points1->
x[vert1 + 1] == Points2->
x[vert2 + 1] &&
1030 Points1->
y[vert1 + 1] == Points2->
y[vert2 + 1]) ||
1031 (Points1->
x[vert1 - 1] == Points2->
x[vert2 + 1] &&
1032 Points1->
y[vert1 - 1] == Points2->
y[vert2 + 1] &&
1033 Points1->
x[vert1 + 1] == Points2->
x[vert2 - 1] &&
1034 Points1->
y[vert1 + 1] == Points2->
y[vert2 - 1])
1037 G_debug(3,
" -> previous/next are not identical");
1043 G_debug(3,
" -> collinear -> remove");
1062 for (i = 0; i < n_cross; i++) {
1063 if (use_cross[i] == 0)
1069 seg = cross[i].segment[current];
1071 G_debug(3,
" duplicate ?: cross = %d seg = %d dist = %f", i,
1072 cross[i].segment[current], cross[i].distance[current]);
1073 if ((cross[i].segment[current] == cross[last].segment[current] &&
1074 cross[i].distance[current] == cross[last].distance[current])
1075 || (cross[i].segment[current] ==
1076 cross[last].segment[current] + 1 &&
1077 cross[i].distance[current] == 0 &&
1078 cross[i].x == cross[last].x &&
1079 cross[i].y == cross[last].y)) {
1080 G_debug(3,
" cross %d identical to last -> removed", i);
1091 G_debug(3,
" alive crosses:");
1092 for (i = 0; i < n_cross; i++) {
1093 if (use_cross[i] == 1) {
1100 if (n_alive_cross > 0) {
1102 use_cross[n_cross] = 1;
1104 cross[n_cross].x = Points->
x[j];
1105 cross[n_cross].y = Points->
y[j];
1106 cross[n_cross].segment[current] = Points->
n_points - 2;
1109 last_x = Points->
x[0];
1110 last_y = Points->
y[0];
1111 last_z = Points->
z[0];
1114 for (i = 0; i <= n_cross; i++) {
1115 seg = cross[i].segment[current];
1116 G_debug(2,
"%d seg = %d dist = %f", i, seg,
1117 cross[i].distance[current]);
1118 if (use_cross[i] == 0) {
1119 G_debug(3,
" removed -> next");
1127 G_debug(2,
" append last vert: %f %f", last_x, last_y);
1130 for (j = last_seg + 1; j <= seg; j++) {
1131 G_debug(2,
" segment j = %d", j);
1133 if ((j == last_seg + 1) && Points->
x[j] == last_x &&
1134 Points->
y[j] == last_y) {
1135 G_debug(2,
" -> skip (identical to last break)");
1140 G_debug(2,
" append first of seg: %f %f", Points->
x[j],
1145 last_x = cross[i].x;
1146 last_y = cross[i].y;
1149 if (Points->
z[last_seg] == Points->
z[last_seg + 1]) {
1150 last_z = Points->
z[last_seg + 1];
1152 else if (last_x == Points->
x[last_seg] &&
1153 last_y == Points->
y[last_seg]) {
1154 last_z = Points->
z[last_seg];
1156 else if (last_x == Points->
x[last_seg + 1] &&
1157 last_y == Points->
y[last_seg + 1]) {
1158 last_z = Points->
z[last_seg + 1];
1161 dist = dist2(Points->
x[last_seg],
1162 Points->
x[last_seg + 1],
1163 Points->
y[last_seg],
1164 Points->
y[last_seg + 1]);
1166 last_z = (Points->
z[last_seg] * sqrt(cross[i].distance[current]) +
1167 Points->
z[last_seg + 1] *
1168 (sqrt(dist) - sqrt(cross[i].distance[current]))) /
1175 G_debug(2,
" append cross / last point: %f %f", cross[i].x,
1180 G_debug(2,
" line is degenerate -> skipped");
1201 static int cross_found;
1204 static int find_cross(
int i,
int j,
int b)
1206 double x1, y1, z1, x2, y2, z2;
1207 double y1min, y1max, y2min, y2max;
1210 y1min = APnts->
y[i];
1211 y1max = APnts->
y[i + 1];
1212 if (APnts->
y[i] > APnts->
y[i + 1]) {
1213 y1min = APnts->
y[i + 1];
1214 y1max = APnts->
y[i];
1217 y2min = BPnts->
y[j];
1218 y2max = BPnts->
y[j + 1];
1219 if (BPnts->
y[j] > BPnts->
y[j + 1]) {
1220 y2min = BPnts->
y[j + 1];
1221 y2max = BPnts->
y[j];
1224 if (y1min > y2max || y1max < y2min)
1230 APnts->
x[i + 1], APnts->
y[i + 1],
1232 BPnts->
x[j], BPnts->
y[j],
1234 BPnts->
x[j + 1], BPnts->
y[j + 1],
1236 &x1, &y1, &z1, &x2, &y2, &z2, 0);
1241 BPnts->
x[j + 1], BPnts->
y[j + 1],
1243 APnts->
x[i], APnts->
y[i],
1245 APnts->
x[i + 1], APnts->
y[i + 1],
1247 &x1, &y1, &z1, &x2, &y2, &z2, 0);
1260 G_warning(
_(
"Error while adding point to array. Out of memory"));
1266 G_warning(
_(
"Error while adding point to array. Out of memory"));
1268 G_warning(
_(
"Error while adding point to array. Out of memory"));
1294 struct boq bo_queue;
1295 struct qitem qi, *found;
1296 struct RB_TREE *bo_ta, *bo_tb;
1299 double xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2, xi, yi;
1315 if (APoints->
x[0] == BPoints->
x[0] && APoints->
y[0] == BPoints->
y[0]) {
1319 &APoints->
y[0],
NULL, 1))
1320 G_warning(
_(
"Error while adding point to array. Out of memory"));
1324 if (APoints->
z[0] == BPoints->
z[0]) {
1327 &APoints->
y[0], &APoints->
z[0],
1329 G_warning(
_(
"Error while adding point to array. Out of memory"));
1346 if (dist <= d_ulp(APoints->
x[0], APoints->
y[0])) {
1350 G_warning(
_(
"Error while adding point to array. Out of memory"));
1363 if (dist <= d_ulp(BPoints->
x[0], BPoints->
y[0])) {
1367 G_warning(
_(
"Error while adding point to array. Out of memory"));
1389 if (abbox.
N > ABox.
N)
1391 if (abbox.
S < ABox.
S)
1393 if (abbox.
E > ABox.
E)
1395 if (abbox.
W < ABox.
W)
1398 abbox.
N += d_ulp(abbox.
N, abbox.
N);
1399 abbox.
S -= d_ulp(abbox.
S, abbox.
S);
1400 abbox.
E += d_ulp(abbox.
E, abbox.
E);
1401 abbox.
W -= d_ulp(abbox.
W, abbox.
W);
1404 if (abbox.
T > ABox.
T)
1406 if (abbox.
B < ABox.
B)
1409 abbox.
T += d_ulp(abbox.
T, abbox.
T);
1410 abbox.
B -= d_ulp(abbox.
B, abbox.
B);
1416 bo_queue.i =
G_malloc(bo_queue.alloc *
sizeof(
struct qitem));
1419 if (!boq_load(&bo_queue, APnts, &abbox, 0, with_z)) {
1425 if (!boq_load(&bo_queue, BPnts, &abbox, 1, with_z)) {
1444 while (boq_drop(&bo_queue, &qi)) {
1452 ret = find_cross(qi.s, found->s, 0);
1459 else if (intersect != 1) {
1463 if (xi == xa1 && yi == ya1) {
1464 if ((xi == xb1 && yi == yb1) ||
1465 (xi == xb2 && yi == yb2)) {
1469 else if (xi == xa2 && yi == ya2) {
1470 if ((xi == xb1 && yi == yb1) ||
1471 (xi == xb2 && yi == yb2)) {
1478 if (intersect == 1) {
1490 ret = find_cross(found->s, qi.s, 1);
1497 else if (intersect != 1) {
1501 if (xi == xa1 && yi == ya1) {
1502 if ((xi == xb1 && yi == yb1) ||
1503 (xi == xb2 && yi == yb2)) {
1507 else if (xi == xa2 && yi == ya2) {
1508 if ((xi == xb1 && yi == yb1) ||
1509 (xi == xb2 && yi == yb2)) {
1516 if (intersect == 1) {
1543 if (intersect == 1) {
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
int dig_line_box(const struct line_pnts *, struct bound_box *)
int n_points
Number of points.
int rbtree_init_trav(struct RB_TRAV *, struct RB_TREE *)
int Vect_line_get_intersections2(struct line_pnts *APoints, struct line_pnts *BPoints, struct line_pnts *IPoints, int with_z)
Get 2 lines intersection 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.
void G_free(void *)
Free allocated memory.
int Vect_line_distance(const struct line_pnts *, double, double, double, int, double *, double *, double *, double *, double *, double *)
Calculate distance of point to line.
struct RB_TREE * rbtree_create(rb_compare_fn *, size_t)
int Vect_segment_intersection(double, double, double, double, double, double, double, double, double, double, double, double, double *, double *, double *, double *, double *, double *, int)
Check for intersect of 2 line segments.
double * x
Array of X coordinates.
Feature geometry info - coordinates.
int rbtree_insert(struct RB_TREE *, void *)
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
void * rbtree_traverse(struct RB_TRAV *)
#define PORT_DOUBLE_MAX
Limits for portable types.
double * y
Array of Y coordinates.
int Vect_line_intersection2(struct line_pnts *APoints, struct line_pnts *BPoints, struct bound_box *pABox, struct bound_box *pBBox, struct line_pnts ***ALines, struct line_pnts ***BLines, int *nalines, int *nblines, int with_z)
Intersect 2 lines.
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_line_check_intersection2(struct line_pnts *APoints, struct line_pnts *BPoints, int with_z)
Check if 2 lines intersect.
int Vect_box_overlap(const struct bound_box *, const struct bound_box *)
Tests for overlap of two boxes.
void rbtree_destroy(struct RB_TREE *)
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
void Vect_reset_line(struct line_pnts *)
Reset line.
int rbtree_remove(struct RB_TREE *, const void *)