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