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);
93 struct line_pnts *BPoints,
int with_z,
int all);
105 static int a_cross = 0;
107 static CROSS *cross =
NULL;
108 static int *use_cross =
NULL;
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) {
140 (CROSS *)
G_realloc((
void *)cross, (a_cross + 101) *
sizeof(CROSS));
142 (
int *)
G_realloc((
void *)use_cross, (a_cross + 101) *
sizeof(
int));
148 " add new cross: aseg/dist = %d/%f bseg/dist = %d/%f, x = %f y = %f",
149 asegment, adistance, bsegment, bdistance,
x, y);
150 cross[n_cross].segment[0] = asegment;
151 cross[n_cross].distance[0] = adistance;
152 cross[n_cross].segment[1] = bsegment;
153 cross[n_cross].distance[1] = bdistance;
154 cross[n_cross].x =
x;
155 cross[n_cross].y = y;
159 static int cmp_cross(
const void *pa,
const void *pb)
161 CROSS *p1 = (CROSS *)pa;
162 CROSS *p2 = (CROSS *)pb;
164 if (p1->segment[current] < p2->segment[current])
166 if (p1->segment[current] > p2->segment[current])
169 if (p1->distance[current] < p2->distance[current])
171 if (p1->distance[current] > p2->distance[current])
176 static double dist2(
double x1,
double y1,
double x2,
double y2)
182 return (dx * dx + dy * dy);
187 static int ident(
double x1,
double y1,
double x2,
double y2,
double thresh)
193 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)
298 APnts->
x[i], APnts->
y[i], APnts->
z[i], APnts->
x[i + 1],
299 APnts->
y[i + 1], APnts->
z[i + 1], BPnts->
x[j], BPnts->
y[j],
300 BPnts->
z[j], BPnts->
x[j + 1], BPnts->
y[j + 1], BPnts->
z[j + 1], &x1,
301 &y1, &z1, &x2, &y2, &z2, 0);
305 BPnts->
x[j], BPnts->
y[j], BPnts->
z[j], BPnts->
x[j + 1],
306 BPnts->
y[j + 1], BPnts->
z[j + 1], APnts->
x[i], APnts->
y[i],
307 APnts->
z[i], APnts->
x[i + 1], APnts->
y[i + 1], APnts->
z[i + 1], &x1,
308 &y1, &z1, &x2, &y2, &z2, 0);
313 G_debug(2,
" -> %d x %d: intersection type = %d", i, j, ret);
315 G_debug(3,
" in %f, %f ", x1, y1);
317 snap_cross(i, &adist, j, &bdist, &x1, &y1);
318 add_cross(i, adist, j, bdist, x1, y1);
320 add_cross(j, bdist, i, adist, x1, y1);
322 else if (ret == 2 || ret == 3 || ret == 4 || ret == 5) {
327 G_debug(3,
" in %f, %f; %f, %f", x1, y1, x2, y2);
328 snap_cross(i, &adist, j, &bdist, &x1, &y1);
329 add_cross(i, adist, j, bdist, x1, y1);
331 add_cross(j, bdist, i, adist, x1, y1);
332 snap_cross(i, &adist, j, &bdist, &x2, &y2);
333 add_cross(i, adist, j, bdist, x2, y2);
335 add_cross(j, bdist, i, adist, x2, y2);
346 #define GET_PARENT(p, c) ((p) = (int)(((c) - 2) / 3 + 1))
347 #define GET_CHILD(c, p) ((c) = (int)(((p) * 3) - 1))
364 static int cmp_q_x(
struct qitem *a,
struct qitem *
b)
366 double x1, y1, z1, x2, y2, z2;
368 x1 = ABPnts[a->l]->
x[a->p];
369 y1 = ABPnts[a->l]->
y[a->p];
370 z1 = ABPnts[a->l]->
z[a->p];
372 x2 = ABPnts[
b->l]->
x[
b->p];
373 y2 = ABPnts[
b->l]->
y[
b->p];
374 z2 = ABPnts[
b->l]->
z[
b->p];
404 static int sift_up(
struct boq *q,
int start)
406 register int parent, child;
417 if (cmp_q_x(&a,
b)) {
419 q->i[child] = q->i[parent];
435 static int boq_add(
struct boq *q,
struct qitem *i)
437 if (q->count + 2 >= q->alloc) {
438 q->alloc = q->count + 100;
439 q->i =
G_realloc(q->i, q->alloc *
sizeof(
struct qitem));
441 q->i[q->count + 1] = *i;
442 sift_up(q, q->count + 1);
450 static int boq_drop(
struct boq *q,
struct qitem *qi)
452 register int child, childr, parent;
471 while (
GET_CHILD(child, parent) <= q->count) {
474 for (childr = child + 1; childr <= q->count && childr < i; childr++) {
483 q->i[parent] = q->i[child];
488 if (parent < q->
count) {
489 q->i[parent] = q->i[q->count];
504 static int cmp_t_y(
const void *aa,
const void *bb)
506 double x1, y1, z1, x2, y2, z2;
507 struct qitem *a = (
struct qitem *)aa;
508 struct qitem *
b = (
struct qitem *)bb;
510 x1 = ABPnts[a->l]->
x[a->p];
511 y1 = ABPnts[a->l]->
y[a->p];
512 z1 = ABPnts[a->l]->
z[a->p];
514 x2 = ABPnts[
b->l]->
x[
b->p];
515 y2 = ABPnts[
b->l]->
y[
b->p];
516 z2 = ABPnts[
b->l]->
z[
b->p];
546 double x1, y1, z1, x2, y2, z2;
563 if (x1 == x2 && y1 == y2 && (!with_z || z1 == z2))
590 box.W -= d_ulp(box.W, box.W);
591 box.S -= d_ulp(box.S, box.S);
592 box.B -= d_ulp(box.B, box.B);
593 box.E += d_ulp(box.E, box.E);
594 box.N += d_ulp(box.N, box.N);
595 box.T += d_ulp(box.T, box.T);
680 struct line_pnts ***BLines,
int *nalines,
681 int *nblines,
int with_z)
683 int i, j, k,
l, nl, last_seg, seg, last;
685 double dist, last_x, last_y, last_z;
688 int seg1, seg2, vert1, vert2;
691 struct qitem qi, *found;
696 if (debug_level == -1) {
700 debug_level = atoi(dstr);
793 if (abbox.
N > ABox.
N)
795 if (abbox.
S < ABox.
S)
797 if (abbox.
E > ABox.
E)
799 if (abbox.
W < ABox.
W)
803 if (abbox.
T > BBox.
T)
805 if (abbox.
B < BBox.
B)
810 abbox.
N += d_ulp(abbox.
N, abbox.
N);
811 abbox.
S -= d_ulp(abbox.
S, abbox.
S);
812 abbox.
E += d_ulp(abbox.
E, abbox.
E);
813 abbox.
W -= d_ulp(abbox.
W, abbox.
W);
815 abbox.
T += d_ulp(abbox.
T, abbox.
T);
816 abbox.
B -= d_ulp(abbox.
B, abbox.
B);
820 G_fatal_error(
"Intersection with points is not yet supported");
827 bo_queue.i =
G_malloc(bo_queue.alloc *
sizeof(
struct qitem));
830 boq_load(&bo_queue, APnts, &abbox, 0, with_z);
834 boq_load(&bo_queue, BPnts, &abbox, 1, with_z);
845 while (boq_drop(&bo_queue, &qi)) {
853 cross_seg(qi.s, found->s, 0);
863 cross_seg(found->s, qi.s, 1);
894 G_debug(2,
"n_cross = %d", n_cross);
904 for (
l = 1;
l < nl;
l++) {
905 for (i = 0; i < n_cross; i++)
912 G_debug(2,
"Clean and create array for line A");
920 G_debug(2,
"Clean and create array for line B");
929 qsort((
void *)cross,
sizeof(
char) * n_cross,
sizeof(CROSS), cmp_cross);
933 if (debug_level > 2) {
934 for (i = 0; i < n_cross; i++) {
937 " cross = %d seg1/dist1 = %d/%f seg2/dist2 = %d/%f x = %f "
939 i, cross[i].segment[current],
940 sqrt(cross[i].distance[current]), cross[i].segment[second],
941 sqrt(cross[i].distance[second]), cross[i].
x, cross[i].y);
946 for (i = 0; i < n_cross; i++) {
947 if (use_cross[i] == 1) {
951 if ((cross[i].segment[current] == 0 &&
952 cross[i].
x == Points1->
x[0] &&
953 cross[i].y == Points1->
y[0]) ||
954 (cross[i].segment[current] == j - 1 &&
955 cross[i].x == Points1->
x[j] &&
956 cross[i].y == Points1->
y[j])) {
958 G_debug(3,
"cross %d deleted (first/last point)", i);
975 for (i = 0; i < n_cross; i++) {
976 if (use_cross[i] == 0)
978 G_debug(3,
" is %d between colinear?", i);
980 seg1 = cross[i].segment[current];
981 seg2 = cross[i].segment[second];
984 if (cross[i].
x == Points1->
x[seg1] &&
985 cross[i].y == Points1->
y[seg1]) {
988 else if (cross[i].
x == Points1->
x[seg1 + 1] &&
989 cross[i].y == Points1->
y[seg1 + 1]) {
993 G_debug(3,
" -> is not vertex on 1. line");
1001 if (cross[i].
x == Points2->
x[seg2] &&
1002 cross[i].y == Points2->
y[seg2]) {
1005 else if (cross[i].
x == Points2->
x[seg2 + 1] &&
1006 cross[i].y == Points2->
y[seg2 + 1]) {
1010 G_debug(3,
" -> is not vertex on 2. line");
1013 G_debug(3,
" seg1/vert1 = %d/%d seg2/ver2 = %d/%d", seg1, vert1,
1017 if (vert2 == 0 || vert2 == Points2->
n_points - 1) {
1018 G_debug(3,
" -> vertex 2 (%d) is first/last", vert2);
1023 if (!((Points1->
x[vert1 - 1] == Points2->
x[vert2 - 1] &&
1024 Points1->
y[vert1 - 1] == Points2->
y[vert2 - 1] &&
1025 Points1->
x[vert1 + 1] == Points2->
x[vert2 + 1] &&
1026 Points1->
y[vert1 + 1] == Points2->
y[vert2 + 1]) ||
1027 (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 G_debug(3,
" -> previous/next are not identical");
1037 G_debug(3,
" -> collinear -> remove");
1058 for (i = 0; i < n_cross; i++) {
1059 if (use_cross[i] == 0)
1065 seg = cross[i].segment[current];
1067 G_debug(3,
" duplicate ?: cross = %d seg = %d dist = %f", i,
1068 cross[i].segment[current], cross[i].distance[current]);
1069 if ((cross[i].segment[current] == cross[last].segment[current] &&
1070 cross[i].distance[current] == cross[last].distance[current]) ||
1071 (cross[i].segment[current] ==
1072 cross[last].segment[current] + 1 &&
1073 cross[i].distance[current] == 0 &&
1074 cross[i].
x == cross[last].
x && cross[i].y == cross[last].y)) {
1075 G_debug(3,
" cross %d identical to last -> removed", i);
1086 G_debug(3,
" alive crosses:");
1087 for (i = 0; i < n_cross; i++) {
1088 if (use_cross[i] == 1) {
1095 if (n_alive_cross > 0) {
1097 use_cross[n_cross] = 1;
1099 cross[n_cross].x = Points->
x[j];
1100 cross[n_cross].y = Points->
y[j];
1101 cross[n_cross].segment[current] = Points->
n_points - 2;
1104 last_x = Points->
x[0];
1105 last_y = Points->
y[0];
1106 last_z = Points->
z[0];
1109 for (i = 0; i <= n_cross; i++) {
1110 seg = cross[i].segment[current];
1111 G_debug(2,
"%d seg = %d dist = %f", i, seg,
1112 cross[i].distance[current]);
1113 if (use_cross[i] == 0) {
1114 G_debug(3,
" removed -> next");
1122 G_debug(2,
" append last vert: %f %f", last_x, last_y);
1125 for (j = last_seg + 1; j <= seg; j++) {
1126 G_debug(2,
" segment j = %d", j);
1128 if ((j == last_seg + 1) && Points->
x[j] == last_x &&
1129 Points->
y[j] == last_y) {
1130 G_debug(2,
" -> skip (identical to last break)");
1135 G_debug(2,
" append first of seg: %f %f", Points->
x[j],
1140 last_x = cross[i].x;
1141 last_y = cross[i].y;
1144 if (Points->
z[last_seg] == Points->
z[last_seg + 1]) {
1145 last_z = Points->
z[last_seg + 1];
1147 else if (last_x == Points->
x[last_seg] &&
1148 last_y == Points->
y[last_seg]) {
1149 last_z = Points->
z[last_seg];
1151 else if (last_x == Points->
x[last_seg + 1] &&
1152 last_y == Points->
y[last_seg + 1]) {
1153 last_z = Points->
z[last_seg + 1];
1156 dist = dist2(Points->
x[last_seg], Points->
x[last_seg + 1],
1157 Points->
y[last_seg], Points->
y[last_seg + 1]);
1159 last_z = (Points->
z[last_seg] *
1160 sqrt(cross[i].distance[current]) +
1161 Points->
z[last_seg + 1] *
1163 sqrt(cross[i].distance[current]))) /
1170 G_debug(2,
" append cross / last point: %f %f", cross[i].
x,
1175 G_debug(2,
" line is degenerate -> skipped");
1199 static int find_cross(
int i,
int j,
int b)
1201 double x1, y1, z1, x2, y2, z2;
1202 double y1min, y1max, y2min, y2max;
1205 y1min = APnts->
y[i];
1206 y1max = APnts->
y[i + 1];
1207 if (APnts->
y[i] > APnts->
y[i + 1]) {
1208 y1min = APnts->
y[i + 1];
1209 y1max = APnts->
y[i];
1212 y2min = BPnts->
y[j];
1213 y2max = BPnts->
y[j + 1];
1214 if (BPnts->
y[j] > BPnts->
y[j + 1]) {
1215 y2min = BPnts->
y[j + 1];
1216 y2max = BPnts->
y[j];
1219 if (y1min > y2max || y1max < y2min)
1224 APnts->
x[i], APnts->
y[i], APnts->
z[i], APnts->
x[i + 1],
1225 APnts->
y[i + 1], APnts->
z[i + 1], BPnts->
x[j], BPnts->
y[j],
1226 BPnts->
z[j], BPnts->
x[j + 1], BPnts->
y[j + 1], BPnts->
z[j + 1], &x1,
1227 &y1, &z1, &x2, &y2, &z2, 0);
1231 BPnts->
x[j], BPnts->
y[j], BPnts->
z[j], BPnts->
x[j + 1],
1232 BPnts->
y[j + 1], BPnts->
z[j + 1], APnts->
x[i], APnts->
y[i],
1233 APnts->
z[i], APnts->
x[i + 1], APnts->
y[i + 1], APnts->
z[i + 1], &x1,
1234 &y1, &z1, &x2, &y2, &z2, 0);
1247 G_warning(
_(
"Error while adding point to array. Out of memory"));
1253 G_warning(
_(
"Error while adding point to array. Out of memory"));
1255 G_warning(
_(
"Error while adding point to array. Out of memory"));
1263 struct line_pnts *BPoints,
int with_z,
int all)
1267 struct boq bo_queue;
1268 struct qitem qi, *found;
1269 struct RB_TREE *bo_ta, *bo_tb;
1272 double xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2, xi, yi;
1289 if (APoints->
x[0] == BPoints->
x[0] && APoints->
y[0] == BPoints->
y[0]) {
1292 APoints->
y[0], APoints->
z[0]))
1294 _(
"Error while adding point to array. Out of memory"));
1298 if (APoints->
z[0] == BPoints->
z[0]) {
1301 APoints->
y[0], APoints->
z[0]))
1302 G_warning(
_(
"Error while adding point to array. Out of "
1319 if (dist <= d_ulp(APoints->
x[0], APoints->
y[0])) {
1321 APoints->
y[0], APoints->
z[0]))
1323 _(
"Error while adding point to array. Out of memory"));
1335 if (dist <= d_ulp(BPoints->
x[0], BPoints->
y[0])) {
1337 BPoints->
y[0], BPoints->
z[0]))
1339 _(
"Error while adding point to array. Out of memory"));
1361 if (abbox.
N > ABox.
N)
1363 if (abbox.
S < ABox.
S)
1365 if (abbox.
E > ABox.
E)
1367 if (abbox.
W < ABox.
W)
1370 abbox.
N += d_ulp(abbox.
N, abbox.
N);
1371 abbox.
S -= d_ulp(abbox.
S, abbox.
S);
1372 abbox.
E += d_ulp(abbox.
E, abbox.
E);
1373 abbox.
W -= d_ulp(abbox.
W, abbox.
W);
1376 if (abbox.
T > ABox.
T)
1378 if (abbox.
B < ABox.
B)
1381 abbox.
T += d_ulp(abbox.
T, abbox.
T);
1382 abbox.
B -= d_ulp(abbox.
B, abbox.
B);
1388 bo_queue.i =
G_malloc(bo_queue.alloc *
sizeof(
struct qitem));
1391 if (!boq_load(&bo_queue, APnts, &abbox, 0, with_z)) {
1397 if (!boq_load(&bo_queue, BPnts, &abbox, 1, with_z)) {
1416 while (boq_drop(&bo_queue, &qi)) {
1424 ret = find_cross(qi.s, found->s, 0);
1431 else if (intersect != 1) {
1435 if (xi == xa1 && yi == ya1) {
1436 if ((xi == xb1 && yi == yb1) ||
1437 (xi == xb2 && yi == yb2)) {
1441 else if (xi == xa2 && yi == ya2) {
1442 if ((xi == xb1 && yi == yb1) ||
1443 (xi == xb2 && yi == yb2)) {
1450 if (intersect == 1) {
1462 ret = find_cross(found->s, qi.s, 1);
1469 else if (intersect != 1) {
1473 if (xi == xa1 && yi == ya1) {
1474 if ((xi == xb1 && yi == yb1) ||
1475 (xi == xb2 && yi == yb2)) {
1479 else if (xi == xa2 && yi == ya2) {
1480 if ((xi == xb1 && yi == yb1) ||
1481 (xi == xb2 && yi == yb2)) {
1488 if (intersect == 1) {
1515 if (!all && intersect == 1) {
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
int G_debug(int, const char *,...) __attribute__((format(printf
const char * G_getenv_nofatal(const char *)
Get environment variable.
int rbtree_insert(struct RB_TREE *, void *)
int rbtree_init_trav(struct RB_TRAV *, struct RB_TREE *)
int rbtree_remove(struct RB_TREE *, const void *)
void * rbtree_traverse(struct RB_TRAV *)
void rbtree_destroy(struct RB_TREE *)
struct RB_TREE * rbtree_create(rb_compare_fn *, size_t)
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
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.
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.
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.
#define PORT_DOUBLE_MAX
Limits for portable types.
int dig_line_degenerate(const struct line_pnts *)
int dig_line_box(const struct line_pnts *, struct bound_box *)
GRASS_INTERPFL_EXPORT struct line_pnts * Pnts
GRASS_INTERPFL_EXPORT int count
int Vect_line_check_intersection2(struct line_pnts *APoints, struct line_pnts *BPoints, int with_z)
Check if 2 lines intersect.
int line_check_intersection2(struct line_pnts *APoints, struct line_pnts *BPoints, int with_z, int all)
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_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.
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.