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);
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
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;
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
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;
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);
327  if (APnts == BPnts)
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);
338  if (APnts == BPnts)
340  snap_cross(i, &adist, j, &bdist, &x2, &y2);
342  if (APnts == BPnts)
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
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 {
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;
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;
650
651  /* event out */
652  qi.e = QEVT_OUT;
653  qi.p = vo;
655
657  }
658
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
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
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