GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
intersect2.c
Go to the documentation of this file.
1 /*!
2  \file lib/vector/Vlib/intersect2.c
3 
4  \brief Vector library - intersection
5 
6  Higher level functions for reading/writing/manipulating vectors.
7 
8  Some parts of code taken from grass50 v.spag/linecros.c
9 
10  Based on the following:
11 
12  <code>
13  (ax2-ax1)r1 - (bx2-bx1)r2 = ax2 - ax1
14  (ay2-ay1)r1 - (by2-by1)r2 = ay2 - ay1
15  </code>
16 
17  Solving for r1 and r2, if r1 and r2 are between 0 and 1, then line
18  segments (ax1,ay1)(ax2,ay2) and (bx1,by1)(bx2,by2) intersect.
19 
20  Intersect 2 line segments.
21 
22  Returns: 0 - do not intersect
23  1 - intersect at one point
24  <pre>
25  \ / \ / \ /
26  \/ \/ \/
27  /\ \
28  / \ \
29  2 - partial overlap ( \/ )
30  ------ a ( distance < threshold )
31  ------ b ( )
32  3 - a contains b ( /\ )
33  ---------- a ----------- a
34  ---- b ----- b
35  4 - b contains a
36  ---- a ----- a
37  ---------- b ----------- b
38  5 - identical
39  ---------- a
40  ---------- b
41  </pre>
42  Intersection points:
43  <pre>
44  return point1 breakes: point2 breaks: distance1 on: distance2 on:
45  0 - - - -
46  1 a,b - a b
47  2 a b a b
48  3 a a a a
49  4 b b b b
50  5 - - - -
51  </pre>
52  Sometimes (often) is important to get the same coordinates for a x
53  b and b x a. To reach this, the segments a,b are 'sorted' at the
54  beginning, so that for the same switched segments, results are
55  identical. (reason is that double values are always rounded because
56  of limited number of decimal places and for different order of
57  coordinates, the results would be different)
58 
59  (C) 2001-2014, 2022 by the GRASS Development Team
60 
61  This program is free software under the GNU General Public License
62  (>=v2). Read the file COPYING that comes with GRASS for details.
63 
64  \author Original author CERL, probably Dave Gerdes or Mike Higgins.
65  \author Update to GRASS 5.7 Radim Blazek.
66  \author Update to GRASS 7 Markus Metz.
67  */
68 
69 #include <stdlib.h>
70 #include <stdio.h>
71 #include <unistd.h>
72 #include <math.h>
73 #include <grass/vector.h>
74 #include <grass/rbtree.h>
75 #include <grass/glocale.h>
76 
77 /* function prototypes */
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);
82 
83 static int debug_level = -1;
84 
85 #if 0
86 static int ident(double x1, double y1, double x2, double y2, double thresh);
87 #endif
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);
92 int line_check_intersection2(struct line_pnts *APoints,
93  struct line_pnts *BPoints, int with_z, int all);
94 
95 typedef struct
96 { /* in arrays 0 - A line , 1 - B line */
97  int segment[2]; /* segment number, start from 0 for first */
98  double distance[2];
99  double x, y, z;
100 } CROSS;
101 
102 /* Current line in arrays is for some functions like cmp() set by: */
103 static int current;
104 static int second; /* line which is not current */
105 
106 static int a_cross = 0;
107 static int n_cross;
108 static CROSS *cross = NULL;
109 static int *use_cross = NULL;
110 
111 /* static double rethresh = 0.000001; */ /* TODO */
112 
113 static double d_ulp(double a, double b)
114 {
115  double fa = fabs(a);
116  double fb = fabs(b);
117  double dmax, result;
118  int exp;
119 
120  dmax = fa;
121  if (dmax < fb)
122  dmax = fb;
123 
124  /* unit in the last place (ULP):
125  * smallest representable difference
126  * shift of the exponent
127  * float: 23, double: 52, middle: 37.5 */
128  result = frexp(dmax, &exp);
129  exp -= 38;
130  result = ldexp(result, exp);
131 
132  return result;
133 }
134 
135 static void add_cross(int asegment, double adistance, int bsegment,
136  double bdistance, double x, double y)
137 {
138  if (n_cross == a_cross) {
139  /* Must be space + 1, used later for last line point, do it better */
140  cross =
141  (CROSS *) G_realloc((void *)cross,
142  (a_cross + 101) * sizeof(CROSS));
143  use_cross =
144  (int *)G_realloc((void *)use_cross,
145  (a_cross + 101) * sizeof(int));
146  a_cross += 100;
147  }
148 
149  G_debug(5,
150  " add new cross: aseg/dist = %d/%f bseg/dist = %d/%f, x = %f y = %f",
151  asegment, adistance, bsegment, bdistance, x, y);
152  cross[n_cross].segment[0] = asegment;
153  cross[n_cross].distance[0] = adistance;
154  cross[n_cross].segment[1] = bsegment;
155  cross[n_cross].distance[1] = bdistance;
156  cross[n_cross].x = x;
157  cross[n_cross].y = y;
158  n_cross++;
159 }
160 
161 static int cmp_cross(const void *pa, const void *pb)
162 {
163  CROSS *p1 = (CROSS *) pa;
164  CROSS *p2 = (CROSS *) pb;
165 
166  if (p1->segment[current] < p2->segment[current])
167  return -1;
168  if (p1->segment[current] > p2->segment[current])
169  return 1;
170  /* the same segment */
171  if (p1->distance[current] < p2->distance[current])
172  return -1;
173  if (p1->distance[current] > p2->distance[current])
174  return 1;
175  return 0;
176 }
177 
178 static double dist2(double x1, double y1, double x2, double y2)
179 {
180  double dx, dy;
181 
182  dx = x2 - x1;
183  dy = y2 - y1;
184  return (dx * dx + dy * dy);
185 }
186 
187 #if 0
188 /* returns 1 if points are identical */
189 static int ident(double x1, double y1, double x2, double y2, double thresh)
190 {
191  double dx, dy;
192 
193  dx = x2 - x1;
194  dy = y2 - y1;
195  if ((dx * dx + dy * dy) <= thresh * thresh)
196  return 1;
197 
198  return 0;
199 }
200 #endif
201 
202 /* shared by Vect_line_intersection, Vect_line_check_intersection, cross_seg, find_cross */
203 static struct line_pnts *APnts, *BPnts, *ABPnts[2], *IPnts;
204 
205 /* Snap breaks to nearest vertices within RE threshold */
206 /* Calculate distances along segments */
207 static int snap_cross(int asegment, double *adistance, int bsegment,
208  double *bdistance, double *xc, double *yc)
209 {
210  int seg;
211  double x, y;
212  double dist, curdist, dthresh;
213 
214  /* 1. of A seg */
215  seg = asegment;
216  curdist = dist2(*xc, *yc, APnts->x[seg], APnts->y[seg]);
217  x = APnts->x[seg];
218  y = APnts->y[seg];
219 
220  *adistance = curdist;
221 
222  /* 2. of A seg */
223  dist = dist2(*xc, *yc, APnts->x[seg + 1], APnts->y[seg + 1]);
224  if (dist < curdist) {
225  curdist = dist;
226  x = APnts->x[seg + 1];
227  y = APnts->y[seg + 1];
228  }
229 
230  /* 1. of B seg */
231  seg = bsegment;
232  dist = dist2(*xc, *yc, BPnts->x[seg], BPnts->y[seg]);
233  *bdistance = dist;
234 
235  if (dist < curdist) {
236  curdist = dist;
237  x = BPnts->x[seg];
238  y = BPnts->y[seg];
239  }
240  /* 2. of B seg */
241  dist = dist2(*xc, *yc, BPnts->x[seg + 1], BPnts->y[seg + 1]);
242  if (dist < curdist) {
243  curdist = dist;
244  x = BPnts->x[seg + 1];
245  y = BPnts->y[seg + 1];
246  }
247 
248  /* the threshold should not be too small, otherwise we get
249  * too many tiny new segments
250  * the threshold should not be too large, otherwise we might
251  * introduce new crossings
252  * the smallest difference representable with
253  * single precision floating point works well with pathological input
254  * regular input is not affected */
255  dthresh = d_ulp(x, y);
256  if (curdist < dthresh * dthresh) { /* was rethresh * rethresh */
257  *xc = x;
258  *yc = y;
259 
260  /* Update distances along segments */
261  seg = asegment;
262  *adistance = dist2(*xc, *yc, APnts->x[seg], APnts->y[seg]);
263  seg = bsegment;
264  *bdistance = dist2(*xc, *yc, BPnts->x[seg], BPnts->y[seg]);
265 
266  return 1;
267  }
268 
269  return 0;
270 }
271 
272 /* break segments */
273 static int cross_seg(int i, int j, int b)
274 {
275  double x1, y1, z1, x2, y2, z2;
276  double y1min, y1max, y2min, y2max;
277  double adist, bdist;
278  int ret;
279 
280  y1min = APnts->y[i];
281  y1max = APnts->y[i + 1];
282  if (APnts->y[i] > APnts->y[i + 1]) {
283  y1min = APnts->y[i + 1];
284  y1max = APnts->y[i];
285  }
286 
287  y2min = BPnts->y[j];
288  y2max = BPnts->y[j + 1];
289  if (BPnts->y[j] > BPnts->y[j + 1]) {
290  y2min = BPnts->y[j + 1];
291  y2max = BPnts->y[j];
292  }
293 
294  if (y1min > y2max || y1max < y2min)
295  return 0;
296 
297  if (b) {
298  ret = Vect_segment_intersection(APnts->x[i], APnts->y[i],
299  APnts->z[i],
300  APnts->x[i + 1], APnts->y[i + 1],
301  APnts->z[i + 1],
302  BPnts->x[j], BPnts->y[j],
303  BPnts->z[j],
304  BPnts->x[j + 1], BPnts->y[j + 1],
305  BPnts->z[j + 1],
306  &x1, &y1, &z1, &x2, &y2, &z2, 0);
307  }
308  else {
309  ret = Vect_segment_intersection(BPnts->x[j], BPnts->y[j],
310  BPnts->z[j],
311  BPnts->x[j + 1], BPnts->y[j + 1],
312  BPnts->z[j + 1],
313  APnts->x[i], APnts->y[i],
314  APnts->z[i],
315  APnts->x[i + 1], APnts->y[i + 1],
316  APnts->z[i + 1],
317  &x1, &y1, &z1, &x2, &y2, &z2, 0);
318  }
319 
320  /* add ALL (including end points and duplicates), clean later */
321  if (ret > 0) {
322  G_debug(2, " -> %d x %d: intersection type = %d", i, j, ret);
323  if (ret == 1) { /* one intersection on segment A */
324  G_debug(3, " in %f, %f ", x1, y1);
325  /* snap intersection only once */
326  snap_cross(i, &adist, j, &bdist, &x1, &y1);
327  add_cross(i, adist, j, bdist, x1, y1);
328  if (APnts == BPnts)
329  add_cross(j, bdist, i, adist, x1, y1);
330  }
331  else if (ret == 2 || ret == 3 || ret == 4 || ret == 5) {
332  /* partial overlap; a broken in one, b broken in one
333  * or a contains b; a is broken in 2 points (but 1 may be end)
334  * or b contains a; b is broken in 2 points (but 1 may be end)
335  * or identical */
336  G_debug(3, " in %f, %f; %f, %f", x1, y1, x2, y2);
337  snap_cross(i, &adist, j, &bdist, &x1, &y1);
338  add_cross(i, adist, j, bdist, x1, y1);
339  if (APnts == BPnts)
340  add_cross(j, bdist, i, adist, x1, y1);
341  snap_cross(i, &adist, j, &bdist, &x2, &y2);
342  add_cross(i, adist, j, bdist, x2, y2);
343  if (APnts == BPnts)
344  add_cross(j, bdist, i, adist, x2, y2);
345  }
346  }
347  return 1; /* keep going */
348 }
349 
350 /* event queue for Bentley-Ottmann */
351 #define QEVT_IN 1
352 #define QEVT_OUT 2
353 #define QEVT_CRS 3
354 
355 #define GET_PARENT(p, c) ((p) = (int) (((c) - 2) / 3 + 1))
356 #define GET_CHILD(c, p) ((c) = (int) (((p) * 3) - 1))
357 
358 struct qitem
359 {
360  int l; /* line 0 - A line , 1 - B line */
361  int s; /* segment index */
362  int p; /* point index */
363  int e; /* event type */
364 };
365 
366 struct boq
367 {
368  int count;
369  int alloc;
370  struct qitem *i;
371 };
372 
373 /* compare two queue points */
374 /* return 1 if a < b else 0 */
375 static int cmp_q_x(struct qitem *a, struct qitem *b)
376 {
377  double x1, y1, z1, x2, y2, z2;
378 
379  x1 = ABPnts[a->l]->x[a->p];
380  y1 = ABPnts[a->l]->y[a->p];
381  z1 = ABPnts[a->l]->z[a->p];
382 
383  x2 = ABPnts[b->l]->x[b->p];
384  y2 = ABPnts[b->l]->y[b->p];
385  z2 = ABPnts[b->l]->z[b->p];
386 
387  if (x1 < x2)
388  return 1;
389  if (x1 > x2)
390  return 0;
391 
392  if (y1 < y2)
393  return 1;
394  if (y1 > y2)
395  return 0;
396 
397  if (z1 < z2)
398  return 1;
399  if (z1 > z2)
400  return 0;
401 
402  if (a->e < b->e)
403  return 1;
404 
405  if (a->l < b->l)
406  return 1;
407 
408  if (a->s < b->s)
409  return 1;
410 
411  return 0;
412 }
413 
414 /* sift up routine for min heap */
415 static int sift_up(struct boq *q, int start)
416 {
417  register int parent, child;
418  struct qitem a, *b;
419 
420  child = start;
421  a = q->i[child];
422 
423  while (child > 1) {
424  GET_PARENT(parent, child);
425 
426  b = &q->i[parent];
427  /* child smaller */
428  if (cmp_q_x(&a, b)) {
429  /* push parent point down */
430  q->i[child] = q->i[parent];
431  child = parent;
432  }
433  else
434  /* no more sifting up, found new slot for child */
435  break;
436  }
437 
438  /* put point in new slot */
439  if (child < start) {
440  q->i[child] = a;
441  }
442 
443  return child;
444 }
445 
446 static int boq_add(struct boq *q, struct qitem *i)
447 {
448  if (q->count + 2 >= q->alloc) {
449  q->alloc = q->count + 100;
450  q->i = G_realloc(q->i, q->alloc * sizeof(struct qitem));
451  }
452  q->i[q->count + 1] = *i;
453  sift_up(q, q->count + 1);
454 
455  q->count++;
456 
457  return 1;
458 }
459 
460 /* drop point routine for min heap */
461 static int boq_drop(struct boq *q, struct qitem *qi)
462 {
463  register int child, childr, parent;
464  register int i;
465  struct qitem *a, *b;
466 
467  if (q->count == 0)
468  return 0;
469 
470  *qi = q->i[1];
471 
472  if (q->count == 1) {
473  q->count = 0;
474  return 1;
475  }
476 
477  /* start with root */
478  parent = 1;
479 
480  /* sift down: move hole back towards bottom of heap */
481 
482  while (GET_CHILD(child, parent) <= q->count) {
483  a = &q->i[child];
484  i = child + 3;
485  for (childr = child + 1; childr <= q->count && childr < i; childr++) {
486  b = &q->i[childr];
487  if (cmp_q_x(b, a)) {
488  child = childr;
489  a = &q->i[child];
490  }
491  }
492 
493  /* move hole down */
494  q->i[parent] = q->i[child];
495  parent = child;
496  }
497 
498  /* hole is in lowest layer, move to heap end */
499  if (parent < q->count) {
500  q->i[parent] = q->i[q->count];
501 
502  /* sift up last swapped point, only necessary if hole moved to heap end */
503  sift_up(q, parent);
504  }
505 
506  /* the actual drop */
507  q->count--;
508 
509  return 1;
510 }
511 
512 /* compare two tree points */
513 /* return -1 if a < b, 1 if a > b, 0 if a == b */
514 static int cmp_t_y(const void *aa, const void *bb)
515 {
516  double x1, y1, z1, x2, y2, z2;
517  struct qitem *a = (struct qitem *) aa;
518  struct qitem *b = (struct qitem *) bb;
519 
520  x1 = ABPnts[a->l]->x[a->p];
521  y1 = ABPnts[a->l]->y[a->p];
522  z1 = ABPnts[a->l]->z[a->p];
523 
524  x2 = ABPnts[b->l]->x[b->p];
525  y2 = ABPnts[b->l]->y[b->p];
526  z2 = ABPnts[b->l]->z[b->p];
527 
528  if (y1 < y2)
529  return -1;
530  if (y1 > y2)
531  return 1;
532 
533  if (x1 < x2)
534  return -1;
535  if (x1 > x2)
536  return 1;
537 
538  if (z1 < z2)
539  return -1;
540  if (z1 > z2)
541  return 1;
542 
543  if (a->s < b->s)
544  return -1;
545  if (a->s > b->s)
546  return 1;
547 
548  return 0;
549 }
550 
551 static int boq_load(struct boq *q, struct line_pnts *Pnts,
552  struct bound_box *abbox, int l, int with_z)
553 {
554  int i, loaded;
555  int vi, vo;
556  double x1, y1, z1, x2, y2, z2;
557  struct bound_box box;
558  struct qitem qi;
559 
560  /* load Pnts to queue */
561  qi.l = l;
562  loaded = 0;
563 
564  for (i = 0; i < Pnts->n_points - 1; i++) {
565  x1 = Pnts->x[i];
566  y1 = Pnts->y[i];
567  z1 = Pnts->z[i];
568 
569  x2 = Pnts->x[i + 1];
570  y2 = Pnts->y[i + 1];
571  z2 = Pnts->z[i + 1];
572 
573  if (x1 == x2 && y1 == y2 && (!with_z || z1 == z2))
574  continue;
575 
576  if (x1 < x2) {
577  box.W = x1;
578  box.E = x2;
579  }
580  else {
581  box.E = x1;
582  box.W = x2;
583  }
584  if (y1 < y2) {
585  box.S = y1;
586  box.N = y2;
587  }
588  else {
589  box.N = y1;
590  box.S = y2;
591  }
592  if (z1 < z2) {
593  box.B = z1;
594  box.T = z2;
595  }
596  else {
597  box.T = z1;
598  box.B = z2;
599  }
600  box.W -= d_ulp(box.W, box.W);
601  box.S -= d_ulp(box.S, box.S);
602  box.B -= d_ulp(box.B, box.B);
603  box.E += d_ulp(box.E, box.E);
604  box.N += d_ulp(box.N, box.N);
605  box.T += d_ulp(box.T, box.T);
606 
607  if (!Vect_box_overlap(abbox, &box))
608  continue;
609 
610  vi = i;
611  vo = i + 1;
612 
613  if (x1 < x2) {
614  vi = i;
615  vo = i + 1;
616  }
617  else if (x1 > x2) {
618  vi = i + 1;
619  vo = i;
620  }
621  else {
622  if (y1 < y2) {
623  vi = i;
624  vo = i + 1;
625  }
626  else if (y1 > y2) {
627  vi = i + 1;
628  vo = i;
629  }
630  else {
631  if (z1 < z2) {
632  vi = i;
633  vo = i + 1;
634  }
635  else if (z1 > z2) {
636  vi = i + 1;
637  vo = i;
638  }
639  else {
640  G_fatal_error("Identical points");
641  }
642  }
643  }
644 
645  qi.s = i;
646 
647  /* event in */
648  qi.e = QEVT_IN;
649  qi.p = vi;
650  boq_add(q, &qi);
651 
652  /* event out */
653  qi.e = QEVT_OUT;
654  qi.p = vo;
655  boq_add(q, &qi);
656 
657  loaded += 2;
658  }
659 
660  return loaded;
661 }
662 
663 /*!
664  * \brief Intersect 2 lines.
665  *
666  * Creates array of new lines created from original A line, by
667  * intersection with B line. Points (Points->n_points == 1) are not
668  * supported. If B line is NULL, A line is intersected with itself.
669  *
670  * simplified Bentley–Ottmann Algorithm:
671  * similar to Vect_line_intersection(), but faster
672  * additionally, self-intersections of a line are handled more efficiently
673  *
674  * \param APoints first input line
675  * \param BPoints second input line or NULL
676  * \param[out] ALines array of new lines created from original A line
677  * \param[out] BLines array of new lines created from original B line
678  * \param[out] nalines number of new lines (ALines)
679  * \param[out] nblines number of new lines (BLines)
680  * \param with_z 3D, not supported!
681  *
682  * \return 0 no intersection
683  * \return 1 intersection found
684  */
685 int
687  struct line_pnts *BPoints,
688  struct bound_box *pABox,
689  struct bound_box *pBBox,
690  struct line_pnts ***ALines,
691  struct line_pnts ***BLines,
692  int *nalines, int *nblines, int with_z)
693 {
694  int i, j, k, l, nl, last_seg, seg, last;
695  int n_alive_cross;
696  double dist, last_x, last_y, last_z;
697  struct line_pnts **XLines, *Points;
698  struct line_pnts *Points1, *Points2; /* first, second points */
699  int seg1, seg2, vert1, vert2;
700  struct bound_box ABox, BBox, abbox;
701  struct boq bo_queue;
702  struct qitem qi, *found;
703  struct RB_TREE *bo_ta, *bo_tb;
704  struct RB_TRAV bo_t_trav;
705  int same = 0;
706 
707  if (debug_level == -1) {
708  const char *dstr = G_getenv_nofatal("DEBUG");
709 
710  if (dstr != NULL)
711  debug_level = atoi(dstr);
712  else
713  debug_level = 0;
714  }
715 
716  n_cross = 0;
717  APnts = APoints;
718  BPnts = BPoints;
719 
720  same = 0;
721  if (!BPoints) {
722  BPnts = APoints;
723  same = 1;
724  }
725 
726  ABPnts[0] = APnts;
727  ABPnts[1] = BPnts;
728 
729  *nalines = 0;
730  *nblines = 0;
731 
732  /* RE (representation error).
733  * RE thresh above is nonsense of course, the RE threshold should be based on
734  * number of significant digits for double (IEEE-754) which is 15 or 16 and exponent.
735  * The number above is in fact not required threshold, and will not work
736  * for example: equator length is 40.075,695 km (8 digits), units are m (+3)
737  * and we want precision in mm (+ 3) = 14 -> minimum rethresh may be around 0.001
738  * ?Maybe all nonsense?
739  * Use rounding error of the unit in the last place ?
740  * max of fabs(x), fabs(y)
741  * rethresh = pow(2, log2(max) - 53) */
742 
743  /* Warning: This function is also used to intersect the line by itself i.e. APoints and
744  * BPoints are identical. I am not sure if it is clever, but it seems to work, but
745  * we have to keep this in mind and handle some special cases (maybe) */
746 
747  /* TODO: 3D, RE threshold, GV_POINTS (line x point) */
748 
749  /* Take each segment from A and intersect by each segment from B.
750  *
751  * All intersections are found first and saved to array, then sorted by a distance along the line,
752  * and then the line is split to pieces.
753  *
754  * Note: If segments are collinear, check if previous/next segments are also collinear,
755  * in that case do not break:
756  * +----------+
757  * +----+-----+ etc.
758  * doesn't need to be broken
759  *
760  * Note: If 2 adjacent segments of line B have common vertex exactly (or within thresh) on line A,
761  * intersection points for these B segments may differ due to RE:
762  * ------------ a ----+--+---- ----+--+----
763  * /\ => / \ or maybe \/
764  * b0 / \ b1 / \ even: /\
765  *
766  * -> solution: snap all breaks to nearest vertices first within RE threshold
767  *
768  * Question: Snap all breaks to each other within RE threshold?
769  *
770  * Note: If a break is snapped to end point or two breaks are snapped to the same vertex
771  * resulting new line is degenerated => before line is added to array, it must be checked
772  * if line is not degenerated
773  *
774  * Note: to snap to vertices is important for cases where line A is broken by B and C line
775  * at the same point:
776  * \ / b no snap \ /
777  * \/ could ----+--+----
778  * ------ a result
779  * /\ in ?: /\
780  * / \ c / \
781  *
782  * Note: once we snap breaks to vertices, we have to do that for both lines A and B in the same way
783  * and because we cannot be sure that A childrens will not change a bit by break(s)
784  * we have to break both A and B at once i.e. in one Vect_line_intersection () call.
785  */
786 
787  /* don't modify original bboxes: make a copy of the bboxes */
788  ABox = *pABox;
789  BBox = *pBBox;
790  if (!with_z) {
791  ABox.T = BBox.T = PORT_DOUBLE_MAX;
792  ABox.B = BBox.B = -PORT_DOUBLE_MAX;
793  }
794 
795  if (!same && !Vect_box_overlap(&ABox, &BBox)) {
796  return 0;
797  }
798 
799  /* overlap box of A line and B line */
800  abbox = BBox;
801  if (!same) {
802  if (abbox.N > ABox.N)
803  abbox.N = ABox.N;
804  if (abbox.S < ABox.S)
805  abbox.S = ABox.S;
806  if (abbox.E > ABox.E)
807  abbox.E = ABox.E;
808  if (abbox.W < ABox.W)
809  abbox.W = ABox.W;
810 
811  if (with_z) {
812  if (abbox.T > BBox.T)
813  abbox.T = BBox.T;
814  if (abbox.B < BBox.B)
815  abbox.B = BBox.B;
816  }
817  }
818 
819  abbox.N += d_ulp(abbox.N, abbox.N);
820  abbox.S -= d_ulp(abbox.S, abbox.S);
821  abbox.E += d_ulp(abbox.E, abbox.E);
822  abbox.W -= d_ulp(abbox.W, abbox.W);
823  if (with_z) {
824  abbox.T += d_ulp(abbox.T, abbox.T);
825  abbox.B -= d_ulp(abbox.B, abbox.B);
826  }
827 
828  if (APnts->n_points < 2 || BPnts->n_points < 2) {
829  G_fatal_error("Intersection with points is not yet supported");
830  return 0;
831  }
832 
833  /* initialize queue */
834  bo_queue.count = 0;
835  bo_queue.alloc = 2 * (APnts->n_points + BPnts->n_points);
836  bo_queue.i = G_malloc(bo_queue.alloc * sizeof(struct qitem));
837 
838  /* load APnts to queue */
839  boq_load(&bo_queue, APnts, &abbox, 0, with_z);
840 
841  if (!same) {
842  /* load BPnts to queue */
843  boq_load(&bo_queue, BPnts, &abbox, 1, with_z);
844  }
845 
846  /* initialize search tree */
847  bo_ta = rbtree_create(cmp_t_y, sizeof(struct qitem));
848  if (same)
849  bo_tb = bo_ta;
850  else
851  bo_tb = rbtree_create(cmp_t_y, sizeof(struct qitem));
852 
853  /* find intersections */
854  while (boq_drop(&bo_queue, &qi)) {
855 
856  if (qi.e == QEVT_IN) {
857  /* not the original Bentley-Ottmann algorithm */
858  if (qi.l == 0) {
859  /* test for intersection of s with all segments in T */
860  rbtree_init_trav(&bo_t_trav, bo_tb);
861  while ((found = rbtree_traverse(&bo_t_trav))) {
862  cross_seg(qi.s, found->s, 0);
863  }
864 
865  /* insert s in T */
866  rbtree_insert(bo_ta, &qi);
867  }
868  else {
869  /* test for intersection of s with all segments in T */
870  rbtree_init_trav(&bo_t_trav, bo_ta);
871  while ((found = rbtree_traverse(&bo_t_trav))) {
872  cross_seg(found->s, qi.s, 1);
873  }
874 
875  /* insert s in T */
876  rbtree_insert(bo_tb, &qi);
877  }
878  }
879  else if (qi.e == QEVT_OUT) {
880  /* remove from T */
881 
882  /* adjust p */
883  if (qi.p == qi.s)
884  qi.p++;
885  else
886  qi.p--;
887 
888  if (qi.l == 0) {
889  if (!rbtree_remove(bo_ta, &qi))
890  G_fatal_error("RB tree error!");
891  }
892  else {
893  if (!rbtree_remove(bo_tb, &qi))
894  G_fatal_error("RB tree error!");
895  }
896  }
897  }
898  G_free(bo_queue.i);
899  rbtree_destroy(bo_ta);
900  if (!same)
901  rbtree_destroy(bo_tb);
902 
903  G_debug(2, "n_cross = %d", n_cross);
904  /* Lines do not cross each other */
905  if (n_cross == 0) {
906  return 0;
907  }
908 
909  /* l = 1 ~ line A, l = 2 ~ line B */
910  nl = 3;
911  if (same)
912  nl = 2;
913  for (l = 1; l < nl; l++) {
914  for (i = 0; i < n_cross; i++)
915  use_cross[i] = 1;
916 
917  /* Create array of lines */
918  XLines = G_malloc((n_cross + 1) * sizeof(struct line_pnts *));
919 
920  if (l == 1) {
921  G_debug(2, "Clean and create array for line A");
922  Points = APnts;
923  Points1 = APnts;
924  Points2 = BPnts;
925  current = 0;
926  second = 1;
927  }
928  else {
929  G_debug(2, "Clean and create array for line B");
930  Points = BPnts;
931  Points1 = BPnts;
932  Points2 = APnts;
933  current = 1;
934  second = 0;
935  }
936 
937  /* Sort points along lines */
938  qsort((void *)cross, sizeof(char) * n_cross, sizeof(CROSS),
939  cmp_cross);
940 
941  /* Print all (raw) breaks */
942  /* avoid loop when not debugging */
943  if (debug_level > 2) {
944  for (i = 0; i < n_cross; i++) {
945  G_debug(3,
946  " cross = %d seg1/dist1 = %d/%f seg2/dist2 = %d/%f x = %f y = %f",
947  i, cross[i].segment[current],
948  sqrt(cross[i].distance[current]),
949  cross[i].segment[second], sqrt(cross[i].distance[second]),
950  cross[i].x, cross[i].y);
951  }
952  }
953 
954  /* Remove breaks on first/last line vertices */
955  for (i = 0; i < n_cross; i++) {
956  if (use_cross[i] == 1) {
957  j = Points1->n_points - 1;
958 
959  /* Note: */
960  if ((cross[i].segment[current] == 0 &&
961  cross[i].x == Points1->x[0] &&
962  cross[i].y == Points1->y[0]) ||
963  (cross[i].segment[current] == j - 1 &&
964  cross[i].x == Points1->x[j] &&
965  cross[i].y == Points1->y[j])) {
966  use_cross[i] = 0; /* first/last */
967  G_debug(3, "cross %d deleted (first/last point)", i);
968  }
969  }
970  }
971 
972  /* Remove breaks with collinear previous and next segments on 1 and 2 */
973  /* Note: breaks with collinear previous and nex must be remove duplicates,
974  * otherwise some cross may be lost. Example (+ is vertex):
975  * B first cross intersections: A/B segment:
976  * | 0/0, 0/1, 1/0, 1/1 - collinear previous and next
977  * AB -----+----+--- A 0/4, 0/5, 1/4, 1/5 - OK
978  * \___|
979  * B
980  * This should not inluence that break is always on first segment, see below (I hope)
981  */
982  /* TODO: this doesn't find identical with breaks on revious/next */
983  for (i = 0; i < n_cross; i++) {
984  if (use_cross[i] == 0)
985  continue;
986  G_debug(3, " is %d between colinear?", i);
987 
988  seg1 = cross[i].segment[current];
989  seg2 = cross[i].segment[second];
990 
991  /* Is it vertex on 1, which? */
992  if (cross[i].x == Points1->x[seg1] &&
993  cross[i].y == Points1->y[seg1]) {
994  vert1 = seg1;
995  }
996  else if (cross[i].x == Points1->x[seg1 + 1] &&
997  cross[i].y == Points1->y[seg1 + 1]) {
998  vert1 = seg1 + 1;
999  }
1000  else {
1001  G_debug(3, " -> is not vertex on 1. line");
1002  continue;
1003  }
1004 
1005  /* Is it vertex on 2, which? */
1006  /* For 1. line it is easy, because breaks on vertex are always at end vertex
1007  * for 2. line we need to find which vertex is on break if any (vert2 starts from 0) */
1008  if (cross[i].x == Points2->x[seg2] &&
1009  cross[i].y == Points2->y[seg2]) {
1010  vert2 = seg2;
1011  }
1012  else if (cross[i].x == Points2->x[seg2 + 1] &&
1013  cross[i].y == Points2->y[seg2 + 1]) {
1014  vert2 = seg2 + 1;
1015  }
1016  else {
1017  G_debug(3, " -> is not vertex on 2. line");
1018  continue;
1019  }
1020  G_debug(3, " seg1/vert1 = %d/%d seg2/ver2 = %d/%d", seg1,
1021  vert1, seg2, vert2);
1022 
1023  /* Check if the second vertex is not first/last */
1024  if (vert2 == 0 || vert2 == Points2->n_points - 1) {
1025  G_debug(3, " -> vertex 2 (%d) is first/last", vert2);
1026  continue;
1027  }
1028 
1029  /* Are there first vertices of this segment identical */
1030  if (!((Points1->x[vert1 - 1] == Points2->x[vert2 - 1] &&
1031  Points1->y[vert1 - 1] == Points2->y[vert2 - 1] &&
1032  Points1->x[vert1 + 1] == Points2->x[vert2 + 1] &&
1033  Points1->y[vert1 + 1] == Points2->y[vert2 + 1]) ||
1034  (Points1->x[vert1 - 1] == Points2->x[vert2 + 1] &&
1035  Points1->y[vert1 - 1] == Points2->y[vert2 + 1] &&
1036  Points1->x[vert1 + 1] == Points2->x[vert2 - 1] &&
1037  Points1->y[vert1 + 1] == Points2->y[vert2 - 1])
1038  )
1039  ) {
1040  G_debug(3, " -> previous/next are not identical");
1041  continue;
1042  }
1043 
1044  use_cross[i] = 0;
1045 
1046  G_debug(3, " -> collinear -> remove");
1047  }
1048 
1049  /* Remove duplicates, i.e. merge all identical breaks to one.
1050  * We must be careful because two points with identical coordinates may be distant if measured along
1051  * the line:
1052  * | Segments b0 and b1 overlap, b0 runs up, b1 down.
1053  * | Two inersections may be merged for a, because they are identical,
1054  * -----+---- a but cannot be merged for b, because both b0 and b1 must be broken.
1055  * | I.e. Breaks on b have identical coordinates, but there are not identical
1056  * b0 | b1 if measured along line b.
1057  *
1058  * -> Breaks may be merged as identical if lay on the same segment, or on vertex connecting
1059  * 2 adjacent segments the points lay on
1060  *
1061  * Note: if duplicate is on a vertex, the break is removed from next segment =>
1062  * break on vertex is always on first segment of this vertex (used below)
1063  */
1064  last = -1;
1065  for (i = 0; i < n_cross; i++) {
1066  if (use_cross[i] == 0)
1067  continue;
1068  if (last == -1) { /* set first alive */
1069  last = i;
1070  continue;
1071  }
1072  seg = cross[i].segment[current];
1073  /* compare with last */
1074  G_debug(3, " duplicate ?: cross = %d seg = %d dist = %f", i,
1075  cross[i].segment[current], cross[i].distance[current]);
1076  if ((cross[i].segment[current] == cross[last].segment[current] &&
1077  cross[i].distance[current] == cross[last].distance[current])
1078  || (cross[i].segment[current] ==
1079  cross[last].segment[current] + 1 &&
1080  cross[i].distance[current] == 0 &&
1081  cross[i].x == cross[last].x &&
1082  cross[i].y == cross[last].y)) {
1083  G_debug(3, " cross %d identical to last -> removed", i);
1084  use_cross[i] = 0; /* identical */
1085  }
1086  else {
1087  last = i;
1088  }
1089  }
1090 
1091  /* Create array of new lines */
1092  /* Count alive crosses */
1093  n_alive_cross = 0;
1094  G_debug(3, " alive crosses:");
1095  for (i = 0; i < n_cross; i++) {
1096  if (use_cross[i] == 1) {
1097  G_debug(3, " %d", i);
1098  n_alive_cross++;
1099  }
1100  }
1101 
1102  k = 0;
1103  if (n_alive_cross > 0) {
1104  /* Add last line point at the end of cross array (cross alley) */
1105  use_cross[n_cross] = 1;
1106  j = Points->n_points - 1;
1107  cross[n_cross].x = Points->x[j];
1108  cross[n_cross].y = Points->y[j];
1109  cross[n_cross].segment[current] = Points->n_points - 2;
1110 
1111  last_seg = 0;
1112  last_x = Points->x[0];
1113  last_y = Points->y[0];
1114  last_z = Points->z[0];
1115  /* Go through all cross (+last line point) and create for each new line
1116  * starting at last_* and ending at cross (last point) */
1117  for (i = 0; i <= n_cross; i++) { /* i.e. n_cross + 1 new lines */
1118  seg = cross[i].segment[current];
1119  G_debug(2, "%d seg = %d dist = %f", i, seg,
1120  cross[i].distance[current]);
1121  if (use_cross[i] == 0) {
1122  G_debug(3, " removed -> next");
1123  continue;
1124  }
1125 
1126  G_debug(2, " New line:");
1127  XLines[k] = Vect_new_line_struct();
1128  /* add last intersection or first point first */
1129  Vect_append_point(XLines[k], last_x, last_y, last_z);
1130  G_debug(2, " append last vert: %f %f", last_x, last_y);
1131 
1132  /* add first points of segments between last and current seg */
1133  for (j = last_seg + 1; j <= seg; j++) {
1134  G_debug(2, " segment j = %d", j);
1135  /* skipp vertex identical to last break */
1136  if ((j == last_seg + 1) && Points->x[j] == last_x &&
1137  Points->y[j] == last_y) {
1138  G_debug(2, " -> skip (identical to last break)");
1139  continue;
1140  }
1141  Vect_append_point(XLines[k], Points->x[j], Points->y[j],
1142  Points->z[j]);
1143  G_debug(2, " append first of seg: %f %f", Points->x[j],
1144  Points->y[j]);
1145  }
1146 
1147  last_seg = seg;
1148  last_x = cross[i].x;
1149  last_y = cross[i].y;
1150  last_z = 0;
1151  /* calculate last_z */
1152  if (Points->z[last_seg] == Points->z[last_seg + 1]) {
1153  last_z = Points->z[last_seg + 1];
1154  }
1155  else if (last_x == Points->x[last_seg] &&
1156  last_y == Points->y[last_seg]) {
1157  last_z = Points->z[last_seg];
1158  }
1159  else if (last_x == Points->x[last_seg + 1] &&
1160  last_y == Points->y[last_seg + 1]) {
1161  last_z = Points->z[last_seg + 1];
1162  }
1163  else {
1164  dist = dist2(Points->x[last_seg],
1165  Points->x[last_seg + 1],
1166  Points->y[last_seg],
1167  Points->y[last_seg + 1]);
1168  if (with_z) {
1169  last_z = (Points->z[last_seg] * sqrt(cross[i].distance[current]) +
1170  Points->z[last_seg + 1] *
1171  (sqrt(dist) - sqrt(cross[i].distance[current]))) /
1172  sqrt(dist);
1173  }
1174  }
1175 
1176  /* add current cross or end point */
1177  Vect_append_point(XLines[k], cross[i].x, cross[i].y, last_z);
1178  G_debug(2, " append cross / last point: %f %f", cross[i].x,
1179  cross[i].y);
1180 
1181  /* Check if line is degenerate */
1182  if (dig_line_degenerate(XLines[k]) > 0) {
1183  G_debug(2, " line is degenerate -> skipped");
1184  Vect_destroy_line_struct(XLines[k]);
1185  }
1186  else {
1187  k++;
1188  }
1189  }
1190  }
1191  if (l == 1) {
1192  *nalines = k;
1193  *ALines = XLines;
1194  }
1195  else {
1196  *nblines = k;
1197  *BLines = XLines;
1198  }
1199  }
1200 
1201  return 1;
1202 }
1203 
1204 /* static int cross_found; */ /* set by find_cross() */
1205 
1206 /* find segment intersection, used by Vect_line_check_intersection2 */
1207 static int find_cross(int i, int j, int b)
1208 {
1209  double x1, y1, z1, x2, y2, z2;
1210  double y1min, y1max, y2min, y2max;
1211  int ret;
1212 
1213  y1min = APnts->y[i];
1214  y1max = APnts->y[i + 1];
1215  if (APnts->y[i] > APnts->y[i + 1]) {
1216  y1min = APnts->y[i + 1];
1217  y1max = APnts->y[i];
1218  }
1219 
1220  y2min = BPnts->y[j];
1221  y2max = BPnts->y[j + 1];
1222  if (BPnts->y[j] > BPnts->y[j + 1]) {
1223  y2min = BPnts->y[j + 1];
1224  y2max = BPnts->y[j];
1225  }
1226 
1227  if (y1min > y2max || y1max < y2min)
1228  return 0;
1229 
1230  if (b) {
1231  ret = Vect_segment_intersection(APnts->x[i], APnts->y[i],
1232  APnts->z[i],
1233  APnts->x[i + 1], APnts->y[i + 1],
1234  APnts->z[i + 1],
1235  BPnts->x[j], BPnts->y[j],
1236  BPnts->z[j],
1237  BPnts->x[j + 1], BPnts->y[j + 1],
1238  BPnts->z[j + 1],
1239  &x1, &y1, &z1, &x2, &y2, &z2, 0);
1240  }
1241  else {
1242  ret = Vect_segment_intersection(BPnts->x[j], BPnts->y[j],
1243  BPnts->z[j],
1244  BPnts->x[j + 1], BPnts->y[j + 1],
1245  BPnts->z[j + 1],
1246  APnts->x[i], APnts->y[i],
1247  APnts->z[i],
1248  APnts->x[i + 1], APnts->y[i + 1],
1249  APnts->z[i + 1],
1250  &x1, &y1, &z1, &x2, &y2, &z2, 0);
1251  }
1252 
1253  if (!IPnts)
1254  IPnts = Vect_new_line_struct();
1255 
1256  /* add ALL (including end points and duplicates), clean later */
1257  switch (ret) {
1258  case 0:
1259  case 5:
1260  break;
1261  case 1:
1262  if (0 > Vect_append_point(IPnts, x1, y1, z1))
1263  G_warning(_("Error while adding point to array. Out of memory"));
1264  break;
1265  case 2:
1266  case 3:
1267  case 4:
1268  if (0 > Vect_append_point(IPnts, x1, y1, z1))
1269  G_warning(_("Error while adding point to array. Out of memory"));
1270  if (0 > Vect_append_point(IPnts, x2, y2, z2))
1271  G_warning(_("Error while adding point to array. Out of memory"));
1272  break;
1273  }
1274 
1275  return ret;
1276 }
1277 
1278 int
1280  struct line_pnts *BPoints, int with_z, int all)
1281 {
1282  double dist;
1283  struct bound_box ABox, BBox, abbox;
1284  struct boq bo_queue;
1285  struct qitem qi, *found;
1286  struct RB_TREE *bo_ta, *bo_tb;
1287  struct RB_TRAV bo_t_trav;
1288  int ret, intersect;
1289  double xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2, xi, yi;
1290 
1291  APnts = APoints;
1292  BPnts = BPoints;
1293 
1294  ABPnts[0] = APnts;
1295  ABPnts[1] = BPnts;
1296 
1297  /* TODO: 3D, RE (representation error) threshold, GV_POINTS (line x point) */
1298 
1299  if (!IPnts)
1300  IPnts = Vect_new_line_struct();
1301  Vect_reset_line(IPnts);
1302 
1303  /* If one or both are point (Points->n_points == 1) */
1304  if (APoints->n_points == 1 && BPoints->n_points == 1) {
1305  if (APoints->x[0] == BPoints->x[0] && APoints->y[0] == BPoints->y[0]) {
1306  if (!with_z) {
1307  if (all && 0 >
1308  Vect_append_point(IPnts, APoints->x[0],
1309  APoints->y[0], APoints->z[0]))
1310  G_warning(_("Error while adding point to array. Out of memory"));
1311  return 1;
1312  }
1313  else {
1314  if (APoints->z[0] == BPoints->z[0]) {
1315  if (all && 0 >
1316  Vect_append_point(IPnts, APoints->x[0],
1317  APoints->y[0], APoints->z[0]))
1318  G_warning(_("Error while adding point to array. Out of memory"));
1319  return 1;
1320  }
1321  else
1322  return 0;
1323  }
1324  }
1325  else {
1326  return 0;
1327  }
1328  }
1329 
1330  if (APoints->n_points == 1) {
1331  Vect_line_distance(BPoints, APoints->x[0], APoints->y[0],
1332  APoints->z[0], with_z, NULL, NULL, NULL, &dist,
1333  NULL, NULL);
1334 
1335  if (dist <= d_ulp(APoints->x[0], APoints->y[0])) {
1336  if (all && 0 >
1337  Vect_append_point(IPnts, APoints->x[0], APoints->y[0],
1338  APoints->z[0]))
1339  G_warning(_("Error while adding point to array. Out of memory"));
1340  return 1;
1341  }
1342  else {
1343  return 0;
1344  }
1345  }
1346 
1347  if (BPoints->n_points == 1) {
1348  Vect_line_distance(APoints, BPoints->x[0], BPoints->y[0],
1349  BPoints->z[0], with_z, NULL, NULL, NULL, &dist,
1350  NULL, NULL);
1351 
1352  if (dist <= d_ulp(BPoints->x[0], BPoints->y[0])) {
1353  if (all && 0 >
1354  Vect_append_point(IPnts, BPoints->x[0], BPoints->y[0],
1355  BPoints->z[0]))
1356  G_warning(_("Error while adding point to array. Out of memory"));
1357  return 1;
1358  }
1359  else
1360  return 0;
1361  }
1362 
1363  /* Take each segment from A and find if intersects any segment from B. */
1364 
1365  dig_line_box(APoints, &ABox);
1366  dig_line_box(BPoints, &BBox);
1367  if (!with_z) {
1368  ABox.T = BBox.T = PORT_DOUBLE_MAX;
1369  ABox.B = BBox.B = -PORT_DOUBLE_MAX;
1370  }
1371 
1372  if (!Vect_box_overlap(&ABox, &BBox)) {
1373  return 0;
1374  }
1375 
1376  /* overlap box */
1377  abbox = BBox;
1378  if (abbox.N > ABox.N)
1379  abbox.N = ABox.N;
1380  if (abbox.S < ABox.S)
1381  abbox.S = ABox.S;
1382  if (abbox.E > ABox.E)
1383  abbox.E = ABox.E;
1384  if (abbox.W < ABox.W)
1385  abbox.W = ABox.W;
1386 
1387  abbox.N += d_ulp(abbox.N, abbox.N);
1388  abbox.S -= d_ulp(abbox.S, abbox.S);
1389  abbox.E += d_ulp(abbox.E, abbox.E);
1390  abbox.W -= d_ulp(abbox.W, abbox.W);
1391 
1392  if (with_z) {
1393  if (abbox.T > ABox.T)
1394  abbox.T = ABox.T;
1395  if (abbox.B < ABox.B)
1396  abbox.B = ABox.B;
1397 
1398  abbox.T += d_ulp(abbox.T, abbox.T);
1399  abbox.B -= d_ulp(abbox.B, abbox.B);
1400  }
1401 
1402  /* initialize queue */
1403  bo_queue.count = 0;
1404  bo_queue.alloc = 2 * (APnts->n_points + BPnts->n_points);
1405  bo_queue.i = G_malloc(bo_queue.alloc * sizeof(struct qitem));
1406 
1407  /* load APnts to queue */
1408  if (!boq_load(&bo_queue, APnts, &abbox, 0, with_z)) {
1409  G_free(bo_queue.i);
1410  return 0;
1411  }
1412 
1413  /* load BPnts to queue */
1414  if (!boq_load(&bo_queue, BPnts, &abbox, 1, with_z)) {
1415  G_free(bo_queue.i);
1416  return 0;
1417  }
1418 
1419  /* initialize search tree */
1420  bo_ta = rbtree_create(cmp_t_y, sizeof(struct qitem));
1421  bo_tb = rbtree_create(cmp_t_y, sizeof(struct qitem));
1422 
1423  /* find intersection */
1424  xa1 = APnts->x[0];
1425  ya1 = APnts->y[0];
1426  xa2 = APnts->x[APnts->n_points - 1];
1427  ya2 = APnts->y[APnts->n_points - 1];
1428  xb1 = BPnts->x[0];
1429  yb1 = BPnts->y[0];
1430  xb2 = BPnts->x[BPnts->n_points - 1];
1431  yb2 = BPnts->y[BPnts->n_points - 1];
1432  intersect = 0;
1433  while (boq_drop(&bo_queue, &qi)) {
1434 
1435  if (qi.e == QEVT_IN) {
1436  /* not the original Bentley-Ottmann algorithm */
1437  if (qi.l == 0) {
1438  /* test for intersection of s with all segments in T */
1439  rbtree_init_trav(&bo_t_trav, bo_tb);
1440  while ((found = rbtree_traverse(&bo_t_trav))) {
1441  ret = find_cross(qi.s, found->s, 0);
1442 
1443  if (ret > 0) {
1444  if (ret != 1) {
1445  intersect = 1;
1446  }
1447  /* intersect at one point */
1448  else if (intersect != 1) {
1449  intersect = 1;
1450  xi = IPnts->x[IPnts->n_points - 1];
1451  yi = IPnts->y[IPnts->n_points - 1];
1452  if (xi == xa1 && yi == ya1) {
1453  if ((xi == xb1 && yi == yb1) ||
1454  (xi == xb2 && yi == yb2)) {
1455  intersect = 2;
1456  }
1457  }
1458  else if (xi == xa2 && yi == ya2) {
1459  if ((xi == xb1 && yi == yb1) ||
1460  (xi == xb2 && yi == yb2)) {
1461  intersect = 2;
1462  }
1463  }
1464  }
1465  }
1466 
1467  if (intersect == 1) {
1468  break;
1469  }
1470  }
1471 
1472  /* insert s in T */
1473  rbtree_insert(bo_ta, &qi);
1474  }
1475  else {
1476  /* test for intersection of s with all segments in T */
1477  rbtree_init_trav(&bo_t_trav, bo_ta);
1478  while ((found = rbtree_traverse(&bo_t_trav))) {
1479  ret = find_cross(found->s, qi.s, 1);
1480 
1481  if (ret > 0) {
1482  if (ret != 1) {
1483  intersect = 1;
1484  }
1485  /* intersect at one point */
1486  else if (intersect != 1) {
1487  intersect = 1;
1488  xi = IPnts->x[IPnts->n_points - 1];
1489  yi = IPnts->y[IPnts->n_points - 1];
1490  if (xi == xa1 && yi == ya1) {
1491  if ((xi == xb1 && yi == yb1) ||
1492  (xi == xb2 && yi == yb2)) {
1493  intersect = 2;
1494  }
1495  }
1496  else if (xi == xa2 && yi == ya2) {
1497  if ((xi == xb1 && yi == yb1) ||
1498  (xi == xb2 && yi == yb2)) {
1499  intersect = 2;
1500  }
1501  }
1502  }
1503  }
1504 
1505  if (intersect == 1) {
1506  break;
1507  }
1508  }
1509 
1510  /* insert s in T */
1511  rbtree_insert(bo_tb, &qi);
1512  }
1513  }
1514  else if (qi.e == QEVT_OUT) {
1515  /* remove from T */
1516 
1517  /* adjust p */
1518  if (qi.p == qi.s)
1519  qi.p++;
1520  else
1521  qi.p--;
1522 
1523  if (qi.l == 0) {
1524  if (!rbtree_remove(bo_ta, &qi))
1525  G_fatal_error("RB tree error!");
1526  }
1527  else {
1528  if (!rbtree_remove(bo_tb, &qi))
1529  G_fatal_error("RB tree error!");
1530  }
1531  }
1532  if (!all && intersect == 1) {
1533  break;
1534  }
1535  }
1536  G_free(bo_queue.i);
1537  rbtree_destroy(bo_ta);
1538  rbtree_destroy(bo_tb);
1539 
1540  return intersect;
1541 }
1542 
1543 /*!
1544  * \brief Check if 2 lines intersect.
1545  *
1546  * Points (Points->n_points == 1) are also supported.
1547  *
1548  * simplified Bentley–Ottmann Algorithm:
1549  * similar to Vect_line_check_intersection(), but faster
1550  *
1551  * \param APoints first input line
1552  * \param BPoints second input line
1553  * \param with_z 3D, not supported (only if one or both are points)!
1554  *
1555  * \return 0 no intersection
1556  * \return 1 intersection
1557  * \return 2 end points only
1558  */
1559 int
1561  struct line_pnts *BPoints, int with_z)
1562 {
1563  return line_check_intersection2(APoints, BPoints, with_z, 0);
1564 }
1565 
1566 /*!
1567  * \brief Get 2 lines intersection points.
1568  *
1569  * A wrapper around Vect_line_check_intersection2() function.
1570  *
1571  * simplified Bentley–Ottmann Algorithm:
1572  * similar to Vect_line_get_intersections(), but faster
1573  *
1574  * \param APoints first input line
1575  * \param BPoints second input line
1576  * \param[out] IPoints output with intersection points
1577  * \param with_z 3D, not supported (only if one or both are points)!
1578  *
1579  * \return 0 no intersection
1580  * \return 1 intersection found
1581  */
1582 int
1584  struct line_pnts *BPoints,
1585  struct line_pnts *IPoints, int with_z)
1586 {
1587  int ret;
1588 
1589  IPnts = IPoints;
1590  ret = line_check_intersection2(APoints, BPoints, with_z, 1);
1591 
1592  return ret;
1593 }
#define G_malloc(n)
Definition: defs/gis.h:112
int line_check_intersection2(struct line_pnts *APoints, struct line_pnts *BPoints, int with_z, int all)
Definition: intersect2.c:1279
Bounding box.
Definition: dig_structs.h:65
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
#define QEVT_IN
Definition: intersect2.c:351
double W
West.
Definition: dig_structs.h:82
int dig_line_box(const struct line_pnts *, struct bound_box *)
int n_points
Number of points.
Definition: dig_structs.h:1692
int rbtree_init_trav(struct RB_TRAV *, struct RB_TREE *)
Definition: rbtree.c:264
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.
Definition: intersect2.c:1583
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
int count
int Vect_line_distance(const struct line_pnts *, double, double, double, int, double *, double *, double *, double *, double *, double *)
Calculate distance of point to line.
Definition: line.c:651
double E
East.
Definition: dig_structs.h:78
#define QEVT_OUT
Definition: intersect2.c:352
struct RB_TREE * rbtree_create(rb_compare_fn *, size_t)
Definition: rbtree.c:50
Definition: rbtree.h:93
#define NULL
Definition: ccmath.h:32
#define x
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.
Definition: rbtree.h:85
double N
North.
Definition: dig_structs.h:70
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
double l
Definition: r_raster.c:39
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
double b
Definition: r_raster.c:39
int rbtree_insert(struct RB_TREE *, void *)
Definition: rbtree.c:74
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition: line.c:149
double B
Bottom.
Definition: dig_structs.h:90
double T
Top.
Definition: dig_structs.h:86
void * rbtree_traverse(struct RB_TRAV *)
Definition: rbtree.c:281
#define GET_CHILD(c, p)
Definition: intersect2.c:356
#define PORT_DOUBLE_MAX
Limits for portable types.
Definition: dig_defines.h:66
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
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.
Definition: intersect2.c:686
void G_warning(const char *,...) __attribute__((format(printf
double S
South.
Definition: dig_structs.h:74
#define G_realloc(p, n)
Definition: defs/gis.h:114
double * z
Array of Z coordinates.
Definition: dig_structs.h:1688
#define _(str)
Definition: glocale.h:10
const char * G_getenv_nofatal(const char *)
Get environment variable.
Definition: env.c:398
int dig_line_degenerate(const struct line_pnts *)
Definition: angle.c:180
int Vect_line_check_intersection2(struct line_pnts *APoints, struct line_pnts *BPoints, int with_z)
Check if 2 lines intersect.
Definition: intersect2.c:1560
int Vect_box_overlap(const struct bound_box *, const struct bound_box *)
Tests for overlap of two boxes.
struct line_pnts * Pnts
void rbtree_destroy(struct RB_TREE *)
Definition: rbtree.c:520
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
Definition: line.c:77
#define GET_PARENT(p, c)
Definition: intersect2.c:355
int G_debug(int, const char *,...) __attribute__((format(printf
void Vect_reset_line(struct line_pnts *)
Reset line.
Definition: line.c:130
int rbtree_remove(struct RB_TREE *, const void *)
Definition: rbtree.c:154