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