GRASS GIS 8 Programmer's Manual  8.5.0dev(2025)-b35369f6b7
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 pABox
669  * \param pBBox
670  * \param[out] nalines number of new lines (ALines)
671  * \param[out] nblines number of new lines (BLines)
672  * \param with_z 3D, not supported!
673  *
674  * \return 0 no intersection
675  * \return 1 intersection found
676  */
677 int Vect_line_intersection2(struct line_pnts *APoints,
678  struct line_pnts *BPoints, struct bound_box *pABox,
679  struct bound_box *pBBox, struct line_pnts ***ALines,
680  struct line_pnts ***BLines, int *nalines,
681  int *nblines, int with_z)
682 {
683  int i, j, k, l, nl, last_seg, seg, last;
684  int n_alive_cross;
685  double dist, last_x, last_y, last_z;
686  struct line_pnts **XLines, *Points;
687  struct line_pnts *Points1, *Points2; /* first, second points */
688  int seg1, seg2, vert1, vert2;
689  struct bound_box ABox, BBox, abbox;
690  struct boq bo_queue;
691  struct qitem qi, *found;
692  struct RB_TREE *bo_ta, *bo_tb;
693  struct RB_TRAV bo_t_trav;
694  int same = 0;
695 
696  if (debug_level == -1) {
697  const char *dstr = G_getenv_nofatal("DEBUG");
698 
699  if (dstr != NULL)
700  debug_level = atoi(dstr);
701  else
702  debug_level = 0;
703  }
704 
705  n_cross = 0;
706  APnts = APoints;
707  BPnts = BPoints;
708 
709  same = 0;
710  if (!BPoints) {
711  BPnts = APoints;
712  same = 1;
713  }
714 
715  ABPnts[0] = APnts;
716  ABPnts[1] = BPnts;
717 
718  *nalines = 0;
719  *nblines = 0;
720 
721  /* RE (representation error).
722  * RE thresh above is nonsense of course, the RE threshold should be based
723  * on number of significant digits for double (IEEE-754) which is 15 or 16
724  * and exponent. The number above is in fact not required threshold, and
725  * will not work for example: equator length is 40.075,695 km (8 digits),
726  * units are m (+3) and we want precision in mm (+ 3) = 14 -> minimum
727  * rethresh may be around 0.001 ?Maybe all nonsense? Use rounding error of
728  * the unit in the last place ? max of fabs(x), fabs(y) rethresh = pow(2,
729  * log2(max) - 53) */
730 
731  /* Warning: This function is also used to intersect the line by itself i.e.
732  * APoints and BPoints are identical. I am not sure if it is clever, but it
733  * seems to work, but we have to keep this in mind and handle some special
734  * cases (maybe) */
735 
736  /* TODO: 3D, RE threshold, GV_POINTS (line x point) */
737 
738  /* Take each segment from A and intersect by each segment from B.
739  *
740  * All intersections are found first and saved to array, then sorted by a
741  * distance along the line, and then the line is split to pieces.
742  *
743  * Note: If segments are collinear, check if previous/next segments are
744  * also collinear, in that case do not break:
745  * +----------+
746  * +----+-----+ etc.
747  * doesn't need to be broken
748  *
749  * Note: If 2 adjacent segments of line B have common vertex exactly (or
750  * within thresh) on line A, intersection points for these B segments may
751  * differ due to RE:
752  * ------------ a ----+--+---- ----+--+----
753  * /\ => / \ or maybe \/
754  * b0 / \ b1 / \ even: /\
755  *
756  * -> solution: snap all breaks to nearest vertices first within RE
757  * threshold
758  *
759  * Question: Snap all breaks to each other within RE threshold?
760  *
761  * Note: If a break is snapped to end point or two breaks are snapped to
762  * the same vertex resulting new line is degenerated => before line is added
763  * to array, it must be checked if line is not degenerated
764  *
765  * Note: to snap to vertices is important for cases where line A is broken
766  * by B and C line at the same point: \ / b no snap \ /
767  * \/ could ----+--+----
768  * ------ a result
769  * /\ in ?: /\
770  * / \ c / \
771  *
772  * Note: once we snap breaks to vertices, we have to do that for both lines
773  * A and B in the same way and because we cannot be sure that A children
774  * will not change a bit by break(s) we have to break both A and B at once
775  * i.e. in one Vect_line_intersection () call.
776  */
777 
778  /* don't modify original bboxes: make a copy of the bboxes */
779  ABox = *pABox;
780  BBox = *pBBox;
781  if (!with_z) {
782  ABox.T = BBox.T = PORT_DOUBLE_MAX;
783  ABox.B = BBox.B = -PORT_DOUBLE_MAX;
784  }
785 
786  if (!same && !Vect_box_overlap(&ABox, &BBox)) {
787  return 0;
788  }
789 
790  /* overlap box of A line and B line */
791  abbox = BBox;
792  if (!same) {
793  if (abbox.N > ABox.N)
794  abbox.N = ABox.N;
795  if (abbox.S < ABox.S)
796  abbox.S = ABox.S;
797  if (abbox.E > ABox.E)
798  abbox.E = ABox.E;
799  if (abbox.W < ABox.W)
800  abbox.W = ABox.W;
801 
802  if (with_z) {
803  if (abbox.T > BBox.T)
804  abbox.T = BBox.T;
805  if (abbox.B < BBox.B)
806  abbox.B = BBox.B;
807  }
808  }
809 
810  abbox.N += d_ulp(abbox.N, abbox.N);
811  abbox.S -= d_ulp(abbox.S, abbox.S);
812  abbox.E += d_ulp(abbox.E, abbox.E);
813  abbox.W -= d_ulp(abbox.W, abbox.W);
814  if (with_z) {
815  abbox.T += d_ulp(abbox.T, abbox.T);
816  abbox.B -= d_ulp(abbox.B, abbox.B);
817  }
818 
819  if (APnts->n_points < 2 || BPnts->n_points < 2) {
820  G_fatal_error("Intersection with points is not yet supported");
821  return 0;
822  }
823 
824  /* initialize queue */
825  bo_queue.count = 0;
826  bo_queue.alloc = 2 * (APnts->n_points + BPnts->n_points);
827  bo_queue.i = G_malloc(bo_queue.alloc * sizeof(struct qitem));
828 
829  /* load APnts to queue */
830  boq_load(&bo_queue, APnts, &abbox, 0, with_z);
831 
832  if (!same) {
833  /* load BPnts to queue */
834  boq_load(&bo_queue, BPnts, &abbox, 1, with_z);
835  }
836 
837  /* initialize search tree */
838  bo_ta = rbtree_create(cmp_t_y, sizeof(struct qitem));
839  if (same)
840  bo_tb = bo_ta;
841  else
842  bo_tb = rbtree_create(cmp_t_y, sizeof(struct qitem));
843 
844  /* find intersections */
845  while (boq_drop(&bo_queue, &qi)) {
846 
847  if (qi.e == QEVT_IN) {
848  /* not the original Bentley-Ottmann algorithm */
849  if (qi.l == 0) {
850  /* test for intersection of s with all segments in T */
851  rbtree_init_trav(&bo_t_trav, bo_tb);
852  while ((found = rbtree_traverse(&bo_t_trav))) {
853  cross_seg(qi.s, found->s, 0);
854  }
855 
856  /* insert s in T */
857  rbtree_insert(bo_ta, &qi);
858  }
859  else {
860  /* test for intersection of s with all segments in T */
861  rbtree_init_trav(&bo_t_trav, bo_ta);
862  while ((found = rbtree_traverse(&bo_t_trav))) {
863  cross_seg(found->s, qi.s, 1);
864  }
865 
866  /* insert s in T */
867  rbtree_insert(bo_tb, &qi);
868  }
869  }
870  else if (qi.e == QEVT_OUT) {
871  /* remove from T */
872 
873  /* adjust p */
874  if (qi.p == qi.s)
875  qi.p++;
876  else
877  qi.p--;
878 
879  if (qi.l == 0) {
880  if (!rbtree_remove(bo_ta, &qi))
881  G_fatal_error("RB tree error!");
882  }
883  else {
884  if (!rbtree_remove(bo_tb, &qi))
885  G_fatal_error("RB tree error!");
886  }
887  }
888  }
889  G_free(bo_queue.i);
890  rbtree_destroy(bo_ta);
891  if (!same)
892  rbtree_destroy(bo_tb);
893 
894  G_debug(2, "n_cross = %d", n_cross);
895  /* Lines do not cross each other */
896  if (n_cross == 0) {
897  return 0;
898  }
899 
900  /* l = 1 ~ line A, l = 2 ~ line B */
901  nl = 3;
902  if (same)
903  nl = 2;
904  for (l = 1; l < nl; l++) {
905  for (i = 0; i < n_cross; i++)
906  use_cross[i] = 1;
907 
908  /* Create array of lines */
909  XLines = G_malloc((n_cross + 1) * sizeof(struct line_pnts *));
910 
911  if (l == 1) {
912  G_debug(2, "Clean and create array for line A");
913  Points = APnts;
914  Points1 = APnts;
915  Points2 = BPnts;
916  current = 0;
917  second = 1;
918  }
919  else {
920  G_debug(2, "Clean and create array for line B");
921  Points = BPnts;
922  Points1 = BPnts;
923  Points2 = APnts;
924  current = 1;
925  second = 0;
926  }
927 
928  /* Sort points along lines */
929  qsort((void *)cross, sizeof(char) * n_cross, sizeof(CROSS), cmp_cross);
930 
931  /* Print all (raw) breaks */
932  /* avoid loop when not debugging */
933  if (debug_level > 2) {
934  for (i = 0; i < n_cross; i++) {
935  G_debug(
936  3,
937  " cross = %d seg1/dist1 = %d/%f seg2/dist2 = %d/%f x = %f "
938  "y = %f",
939  i, cross[i].segment[current],
940  sqrt(cross[i].distance[current]), cross[i].segment[second],
941  sqrt(cross[i].distance[second]), cross[i].x, cross[i].y);
942  }
943  }
944 
945  /* Remove breaks on first/last line vertices */
946  for (i = 0; i < n_cross; i++) {
947  if (use_cross[i] == 1) {
948  j = Points1->n_points - 1;
949 
950  /* Note: */
951  if ((cross[i].segment[current] == 0 &&
952  cross[i].x == Points1->x[0] &&
953  cross[i].y == Points1->y[0]) ||
954  (cross[i].segment[current] == j - 1 &&
955  cross[i].x == Points1->x[j] &&
956  cross[i].y == Points1->y[j])) {
957  use_cross[i] = 0; /* first/last */
958  G_debug(3, "cross %d deleted (first/last point)", i);
959  }
960  }
961  }
962 
963  /* Remove breaks with collinear previous and next segments on 1 and 2 */
964  /* Note: breaks with collinear previous and nex must be remove
965  * duplicates, otherwise some cross may be lost. Example (+ is vertex):
966  * B first cross intersections: A/B segment:
967  * | 0/0, 0/1, 1/0, 1/1 - collinear previous
968  * and next AB -----+----+--- A 0/4, 0/5, 1/4, 1/5 - OK
969  * \___|
970  * B
971  * This should not influence that break is always on first segment, see
972  * below (I hope)
973  */
974  /* TODO: this doesn't find identical with breaks on revious/next */
975  for (i = 0; i < n_cross; i++) {
976  if (use_cross[i] == 0)
977  continue;
978  G_debug(3, " is %d between colinear?", i);
979 
980  seg1 = cross[i].segment[current];
981  seg2 = cross[i].segment[second];
982 
983  /* Is it vertex on 1, which? */
984  if (cross[i].x == Points1->x[seg1] &&
985  cross[i].y == Points1->y[seg1]) {
986  vert1 = seg1;
987  }
988  else if (cross[i].x == Points1->x[seg1 + 1] &&
989  cross[i].y == Points1->y[seg1 + 1]) {
990  vert1 = seg1 + 1;
991  }
992  else {
993  G_debug(3, " -> is not vertex on 1. line");
994  continue;
995  }
996 
997  /* Is it vertex on 2, which? */
998  /* For 1. line it is easy, because breaks on vertex are always at
999  * end vertex for 2. line we need to find which vertex is on break
1000  * if any (vert2 starts from 0) */
1001  if (cross[i].x == Points2->x[seg2] &&
1002  cross[i].y == Points2->y[seg2]) {
1003  vert2 = seg2;
1004  }
1005  else if (cross[i].x == Points2->x[seg2 + 1] &&
1006  cross[i].y == Points2->y[seg2 + 1]) {
1007  vert2 = seg2 + 1;
1008  }
1009  else {
1010  G_debug(3, " -> is not vertex on 2. line");
1011  continue;
1012  }
1013  G_debug(3, " seg1/vert1 = %d/%d seg2/ver2 = %d/%d", seg1, vert1,
1014  seg2, vert2);
1015 
1016  /* Check if the second vertex is not first/last */
1017  if (vert2 == 0 || vert2 == Points2->n_points - 1) {
1018  G_debug(3, " -> vertex 2 (%d) is first/last", vert2);
1019  continue;
1020  }
1021 
1022  /* Are there first vertices of this segment identical */
1023  if (!((Points1->x[vert1 - 1] == Points2->x[vert2 - 1] &&
1024  Points1->y[vert1 - 1] == Points2->y[vert2 - 1] &&
1025  Points1->x[vert1 + 1] == Points2->x[vert2 + 1] &&
1026  Points1->y[vert1 + 1] == Points2->y[vert2 + 1]) ||
1027  (Points1->x[vert1 - 1] == Points2->x[vert2 + 1] &&
1028  Points1->y[vert1 - 1] == Points2->y[vert2 + 1] &&
1029  Points1->x[vert1 + 1] == Points2->x[vert2 - 1] &&
1030  Points1->y[vert1 + 1] == Points2->y[vert2 - 1]))) {
1031  G_debug(3, " -> previous/next are not identical");
1032  continue;
1033  }
1034 
1035  use_cross[i] = 0;
1036 
1037  G_debug(3, " -> collinear -> remove");
1038  }
1039 
1040  /* Remove duplicates, i.e. merge all identical breaks to one.
1041  * We must be careful because two points with identical coordinates may
1042  * be distant if measured along the line: | Segments b0 and b1
1043  * overlap, b0 runs up, b1 down. | Two inersections may be
1044  * merged for a, because they are identical,
1045  * -----+---- a but cannot be merged for b, because both b0 and b1
1046  * must be broken. | I.e. Breaks on b have identical
1047  * coordinates, but there are not identical b0 | b1 if measured
1048  * along line b.
1049  *
1050  * -> Breaks may be merged as identical if lay on the same segment,
1051  * or on vertex connecting 2 adjacent segments the points lay on
1052  *
1053  * Note: if duplicate is on a vertex, the break is removed from next
1054  * segment => break on vertex is always on first segment of this vertex
1055  * (used below)
1056  */
1057  last = -1;
1058  for (i = 0; i < n_cross; i++) {
1059  if (use_cross[i] == 0)
1060  continue;
1061  if (last == -1) { /* set first alive */
1062  last = i;
1063  continue;
1064  }
1065  seg = cross[i].segment[current];
1066  /* compare with last */
1067  G_debug(3, " duplicate ?: cross = %d seg = %d dist = %f", i,
1068  cross[i].segment[current], cross[i].distance[current]);
1069  if ((cross[i].segment[current] == cross[last].segment[current] &&
1070  cross[i].distance[current] == cross[last].distance[current]) ||
1071  (cross[i].segment[current] ==
1072  cross[last].segment[current] + 1 &&
1073  cross[i].distance[current] == 0 &&
1074  cross[i].x == cross[last].x && cross[i].y == cross[last].y)) {
1075  G_debug(3, " cross %d identical to last -> removed", i);
1076  use_cross[i] = 0; /* identical */
1077  }
1078  else {
1079  last = i;
1080  }
1081  }
1082 
1083  /* Create array of new lines */
1084  /* Count alive crosses */
1085  n_alive_cross = 0;
1086  G_debug(3, " alive crosses:");
1087  for (i = 0; i < n_cross; i++) {
1088  if (use_cross[i] == 1) {
1089  G_debug(3, " %d", i);
1090  n_alive_cross++;
1091  }
1092  }
1093 
1094  k = 0;
1095  if (n_alive_cross > 0) {
1096  /* Add last line point at the end of cross array (cross alley) */
1097  use_cross[n_cross] = 1;
1098  j = Points->n_points - 1;
1099  cross[n_cross].x = Points->x[j];
1100  cross[n_cross].y = Points->y[j];
1101  cross[n_cross].segment[current] = Points->n_points - 2;
1102 
1103  last_seg = 0;
1104  last_x = Points->x[0];
1105  last_y = Points->y[0];
1106  last_z = Points->z[0];
1107  /* Go through all cross (+last line point) and create for each new
1108  * line starting at last_* and ending at cross (last point) */
1109  for (i = 0; i <= n_cross; i++) { /* i.e. n_cross + 1 new lines */
1110  seg = cross[i].segment[current];
1111  G_debug(2, "%d seg = %d dist = %f", i, seg,
1112  cross[i].distance[current]);
1113  if (use_cross[i] == 0) {
1114  G_debug(3, " removed -> next");
1115  continue;
1116  }
1117 
1118  G_debug(2, " New line:");
1119  XLines[k] = Vect_new_line_struct();
1120  /* add last intersection or first point first */
1121  Vect_append_point(XLines[k], last_x, last_y, last_z);
1122  G_debug(2, " append last vert: %f %f", last_x, last_y);
1123 
1124  /* add first points of segments between last and current seg */
1125  for (j = last_seg + 1; j <= seg; j++) {
1126  G_debug(2, " segment j = %d", j);
1127  /* skip vertex identical to last break */
1128  if ((j == last_seg + 1) && Points->x[j] == last_x &&
1129  Points->y[j] == last_y) {
1130  G_debug(2, " -> skip (identical to last break)");
1131  continue;
1132  }
1133  Vect_append_point(XLines[k], Points->x[j], Points->y[j],
1134  Points->z[j]);
1135  G_debug(2, " append first of seg: %f %f", Points->x[j],
1136  Points->y[j]);
1137  }
1138 
1139  last_seg = seg;
1140  last_x = cross[i].x;
1141  last_y = cross[i].y;
1142  last_z = 0;
1143  /* calculate last_z */
1144  if (Points->z[last_seg] == Points->z[last_seg + 1]) {
1145  last_z = Points->z[last_seg + 1];
1146  }
1147  else if (last_x == Points->x[last_seg] &&
1148  last_y == Points->y[last_seg]) {
1149  last_z = Points->z[last_seg];
1150  }
1151  else if (last_x == Points->x[last_seg + 1] &&
1152  last_y == Points->y[last_seg + 1]) {
1153  last_z = Points->z[last_seg + 1];
1154  }
1155  else {
1156  dist = dist2(Points->x[last_seg], Points->x[last_seg + 1],
1157  Points->y[last_seg], Points->y[last_seg + 1]);
1158  if (with_z) {
1159  last_z = (Points->z[last_seg] *
1160  sqrt(cross[i].distance[current]) +
1161  Points->z[last_seg + 1] *
1162  (sqrt(dist) -
1163  sqrt(cross[i].distance[current]))) /
1164  sqrt(dist);
1165  }
1166  }
1167 
1168  /* add current cross or end point */
1169  Vect_append_point(XLines[k], cross[i].x, cross[i].y, last_z);
1170  G_debug(2, " append cross / last point: %f %f", cross[i].x,
1171  cross[i].y);
1172 
1173  /* Check if line is degenerate */
1174  if (dig_line_degenerate(XLines[k]) > 0) {
1175  G_debug(2, " line is degenerate -> skipped");
1176  Vect_destroy_line_struct(XLines[k]);
1177  }
1178  else {
1179  k++;
1180  }
1181  }
1182  }
1183  if (l == 1) {
1184  *nalines = k;
1185  *ALines = XLines;
1186  }
1187  else {
1188  *nblines = k;
1189  *BLines = XLines;
1190  }
1191  }
1192 
1193  return 1;
1194 }
1195 
1196 /* static int cross_found; */ /* set by find_cross() */
1197 
1198 /* find segment intersection, used by Vect_line_check_intersection2 */
1199 static int find_cross(int i, int j, int b)
1200 {
1201  double x1, y1, z1, x2, y2, z2;
1202  double y1min, y1max, y2min, y2max;
1203  int ret;
1204 
1205  y1min = APnts->y[i];
1206  y1max = APnts->y[i + 1];
1207  if (APnts->y[i] > APnts->y[i + 1]) {
1208  y1min = APnts->y[i + 1];
1209  y1max = APnts->y[i];
1210  }
1211 
1212  y2min = BPnts->y[j];
1213  y2max = BPnts->y[j + 1];
1214  if (BPnts->y[j] > BPnts->y[j + 1]) {
1215  y2min = BPnts->y[j + 1];
1216  y2max = BPnts->y[j];
1217  }
1218 
1219  if (y1min > y2max || y1max < y2min)
1220  return 0;
1221 
1222  if (b) {
1224  APnts->x[i], APnts->y[i], APnts->z[i], APnts->x[i + 1],
1225  APnts->y[i + 1], APnts->z[i + 1], BPnts->x[j], BPnts->y[j],
1226  BPnts->z[j], BPnts->x[j + 1], BPnts->y[j + 1], BPnts->z[j + 1], &x1,
1227  &y1, &z1, &x2, &y2, &z2, 0);
1228  }
1229  else {
1231  BPnts->x[j], BPnts->y[j], BPnts->z[j], BPnts->x[j + 1],
1232  BPnts->y[j + 1], BPnts->z[j + 1], APnts->x[i], APnts->y[i],
1233  APnts->z[i], APnts->x[i + 1], APnts->y[i + 1], APnts->z[i + 1], &x1,
1234  &y1, &z1, &x2, &y2, &z2, 0);
1235  }
1236 
1237  if (!IPnts)
1238  IPnts = Vect_new_line_struct();
1239 
1240  /* add ALL (including end points and duplicates), clean later */
1241  switch (ret) {
1242  case 0:
1243  case 5:
1244  break;
1245  case 1:
1246  if (0 > Vect_append_point(IPnts, x1, y1, z1))
1247  G_warning(_("Error while adding point to array. Out of memory"));
1248  break;
1249  case 2:
1250  case 3:
1251  case 4:
1252  if (0 > Vect_append_point(IPnts, x1, y1, z1))
1253  G_warning(_("Error while adding point to array. Out of memory"));
1254  if (0 > Vect_append_point(IPnts, x2, y2, z2))
1255  G_warning(_("Error while adding point to array. Out of memory"));
1256  break;
1257  }
1258 
1259  return ret;
1260 }
1261 
1263  struct line_pnts *BPoints, int with_z, int all)
1264 {
1265  double dist;
1266  struct bound_box ABox, BBox, abbox;
1267  struct boq bo_queue;
1268  struct qitem qi, *found;
1269  struct RB_TREE *bo_ta, *bo_tb;
1270  struct RB_TRAV bo_t_trav;
1271  int ret, intersect;
1272  double xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2, xi, yi;
1273 
1274  APnts = APoints;
1275  BPnts = BPoints;
1276 
1277  ABPnts[0] = APnts;
1278  ABPnts[1] = BPnts;
1279 
1280  /* TODO: 3D, RE (representation error) threshold, GV_POINTS (line x point)
1281  */
1282 
1283  if (!IPnts)
1284  IPnts = Vect_new_line_struct();
1285  Vect_reset_line(IPnts);
1286 
1287  /* If one or both are point (Points->n_points == 1) */
1288  if (APoints->n_points == 1 && BPoints->n_points == 1) {
1289  if (APoints->x[0] == BPoints->x[0] && APoints->y[0] == BPoints->y[0]) {
1290  if (!with_z) {
1291  if (all && 0 > Vect_append_point(IPnts, APoints->x[0],
1292  APoints->y[0], APoints->z[0]))
1293  G_warning(
1294  _("Error while adding point to array. Out of memory"));
1295  return 1;
1296  }
1297  else {
1298  if (APoints->z[0] == BPoints->z[0]) {
1299  if (all &&
1300  0 > Vect_append_point(IPnts, APoints->x[0],
1301  APoints->y[0], APoints->z[0]))
1302  G_warning(_("Error while adding point to array. Out of "
1303  "memory"));
1304  return 1;
1305  }
1306  else
1307  return 0;
1308  }
1309  }
1310  else {
1311  return 0;
1312  }
1313  }
1314 
1315  if (APoints->n_points == 1) {
1316  Vect_line_distance(BPoints, APoints->x[0], APoints->y[0], APoints->z[0],
1317  with_z, NULL, NULL, NULL, &dist, NULL, NULL);
1318 
1319  if (dist <= d_ulp(APoints->x[0], APoints->y[0])) {
1320  if (all && 0 > Vect_append_point(IPnts, APoints->x[0],
1321  APoints->y[0], APoints->z[0]))
1322  G_warning(
1323  _("Error while adding point to array. Out of memory"));
1324  return 1;
1325  }
1326  else {
1327  return 0;
1328  }
1329  }
1330 
1331  if (BPoints->n_points == 1) {
1332  Vect_line_distance(APoints, BPoints->x[0], BPoints->y[0], BPoints->z[0],
1333  with_z, NULL, NULL, NULL, &dist, NULL, NULL);
1334 
1335  if (dist <= d_ulp(BPoints->x[0], BPoints->y[0])) {
1336  if (all && 0 > Vect_append_point(IPnts, BPoints->x[0],
1337  BPoints->y[0], BPoints->z[0]))
1338  G_warning(
1339  _("Error while adding point to array. Out of memory"));
1340  return 1;
1341  }
1342  else
1343  return 0;
1344  }
1345 
1346  /* Take each segment from A and find if intersects any segment from B. */
1347 
1348  dig_line_box(APoints, &ABox);
1349  dig_line_box(BPoints, &BBox);
1350  if (!with_z) {
1351  ABox.T = BBox.T = PORT_DOUBLE_MAX;
1352  ABox.B = BBox.B = -PORT_DOUBLE_MAX;
1353  }
1354 
1355  if (!Vect_box_overlap(&ABox, &BBox)) {
1356  return 0;
1357  }
1358 
1359  /* overlap box */
1360  abbox = BBox;
1361  if (abbox.N > ABox.N)
1362  abbox.N = ABox.N;
1363  if (abbox.S < ABox.S)
1364  abbox.S = ABox.S;
1365  if (abbox.E > ABox.E)
1366  abbox.E = ABox.E;
1367  if (abbox.W < ABox.W)
1368  abbox.W = ABox.W;
1369 
1370  abbox.N += d_ulp(abbox.N, abbox.N);
1371  abbox.S -= d_ulp(abbox.S, abbox.S);
1372  abbox.E += d_ulp(abbox.E, abbox.E);
1373  abbox.W -= d_ulp(abbox.W, abbox.W);
1374 
1375  if (with_z) {
1376  if (abbox.T > ABox.T)
1377  abbox.T = ABox.T;
1378  if (abbox.B < ABox.B)
1379  abbox.B = ABox.B;
1380 
1381  abbox.T += d_ulp(abbox.T, abbox.T);
1382  abbox.B -= d_ulp(abbox.B, abbox.B);
1383  }
1384 
1385  /* initialize queue */
1386  bo_queue.count = 0;
1387  bo_queue.alloc = 2 * (APnts->n_points + BPnts->n_points);
1388  bo_queue.i = G_malloc(bo_queue.alloc * sizeof(struct qitem));
1389 
1390  /* load APnts to queue */
1391  if (!boq_load(&bo_queue, APnts, &abbox, 0, with_z)) {
1392  G_free(bo_queue.i);
1393  return 0;
1394  }
1395 
1396  /* load BPnts to queue */
1397  if (!boq_load(&bo_queue, BPnts, &abbox, 1, with_z)) {
1398  G_free(bo_queue.i);
1399  return 0;
1400  }
1401 
1402  /* initialize search tree */
1403  bo_ta = rbtree_create(cmp_t_y, sizeof(struct qitem));
1404  bo_tb = rbtree_create(cmp_t_y, sizeof(struct qitem));
1405 
1406  /* find intersection */
1407  xa1 = APnts->x[0];
1408  ya1 = APnts->y[0];
1409  xa2 = APnts->x[APnts->n_points - 1];
1410  ya2 = APnts->y[APnts->n_points - 1];
1411  xb1 = BPnts->x[0];
1412  yb1 = BPnts->y[0];
1413  xb2 = BPnts->x[BPnts->n_points - 1];
1414  yb2 = BPnts->y[BPnts->n_points - 1];
1415  intersect = 0;
1416  while (boq_drop(&bo_queue, &qi)) {
1417 
1418  if (qi.e == QEVT_IN) {
1419  /* not the original Bentley-Ottmann algorithm */
1420  if (qi.l == 0) {
1421  /* test for intersection of s with all segments in T */
1422  rbtree_init_trav(&bo_t_trav, bo_tb);
1423  while ((found = rbtree_traverse(&bo_t_trav))) {
1424  ret = find_cross(qi.s, found->s, 0);
1425 
1426  if (ret > 0) {
1427  if (ret != 1) {
1428  intersect = 1;
1429  }
1430  /* intersect at one point */
1431  else if (intersect != 1) {
1432  intersect = 1;
1433  xi = IPnts->x[IPnts->n_points - 1];
1434  yi = IPnts->y[IPnts->n_points - 1];
1435  if (xi == xa1 && yi == ya1) {
1436  if ((xi == xb1 && yi == yb1) ||
1437  (xi == xb2 && yi == yb2)) {
1438  intersect = 2;
1439  }
1440  }
1441  else if (xi == xa2 && yi == ya2) {
1442  if ((xi == xb1 && yi == yb1) ||
1443  (xi == xb2 && yi == yb2)) {
1444  intersect = 2;
1445  }
1446  }
1447  }
1448  }
1449 
1450  if (intersect == 1) {
1451  break;
1452  }
1453  }
1454 
1455  /* insert s in T */
1456  rbtree_insert(bo_ta, &qi);
1457  }
1458  else {
1459  /* test for intersection of s with all segments in T */
1460  rbtree_init_trav(&bo_t_trav, bo_ta);
1461  while ((found = rbtree_traverse(&bo_t_trav))) {
1462  ret = find_cross(found->s, qi.s, 1);
1463 
1464  if (ret > 0) {
1465  if (ret != 1) {
1466  intersect = 1;
1467  }
1468  /* intersect at one point */
1469  else if (intersect != 1) {
1470  intersect = 1;
1471  xi = IPnts->x[IPnts->n_points - 1];
1472  yi = IPnts->y[IPnts->n_points - 1];
1473  if (xi == xa1 && yi == ya1) {
1474  if ((xi == xb1 && yi == yb1) ||
1475  (xi == xb2 && yi == yb2)) {
1476  intersect = 2;
1477  }
1478  }
1479  else if (xi == xa2 && yi == ya2) {
1480  if ((xi == xb1 && yi == yb1) ||
1481  (xi == xb2 && yi == yb2)) {
1482  intersect = 2;
1483  }
1484  }
1485  }
1486  }
1487 
1488  if (intersect == 1) {
1489  break;
1490  }
1491  }
1492 
1493  /* insert s in T */
1494  rbtree_insert(bo_tb, &qi);
1495  }
1496  }
1497  else if (qi.e == QEVT_OUT) {
1498  /* remove from T */
1499 
1500  /* adjust p */
1501  if (qi.p == qi.s)
1502  qi.p++;
1503  else
1504  qi.p--;
1505 
1506  if (qi.l == 0) {
1507  if (!rbtree_remove(bo_ta, &qi))
1508  G_fatal_error("RB tree error!");
1509  }
1510  else {
1511  if (!rbtree_remove(bo_tb, &qi))
1512  G_fatal_error("RB tree error!");
1513  }
1514  }
1515  if (!all && intersect == 1) {
1516  break;
1517  }
1518  }
1519  G_free(bo_queue.i);
1520  rbtree_destroy(bo_ta);
1521  rbtree_destroy(bo_tb);
1522 
1523  return intersect;
1524 }
1525 
1526 /*!
1527  * \brief Check if 2 lines intersect.
1528  *
1529  * Points (Points->n_points == 1) are also supported.
1530  *
1531  * simplified Bentley–Ottmann Algorithm:
1532  * similar to Vect_line_check_intersection(), but faster
1533  *
1534  * \param APoints first input line
1535  * \param BPoints second input line
1536  * \param with_z 3D, not supported (only if one or both are points)!
1537  *
1538  * \return 0 no intersection
1539  * \return 1 intersection
1540  * \return 2 end points only
1541  */
1543  struct line_pnts *BPoints, int with_z)
1544 {
1545  return line_check_intersection2(APoints, BPoints, with_z, 0);
1546 }
1547 
1548 /*!
1549  * \brief Get 2 lines intersection points.
1550  *
1551  * A wrapper around Vect_line_check_intersection2() function.
1552  *
1553  * simplified Bentley–Ottmann Algorithm:
1554  * similar to Vect_line_get_intersections(), but faster
1555  *
1556  * \param APoints first input line
1557  * \param BPoints second input line
1558  * \param[out] IPoints output with intersection points
1559  * \param with_z 3D, not supported (only if one or both are points)!
1560  *
1561  * \return 0 no intersection
1562  * \return 1 intersection found
1563  */
1565  struct line_pnts *BPoints,
1566  struct line_pnts *IPoints, int with_z)
1567 {
1568  int ret;
1569 
1570  IPnts = IPoints;
1571  ret = line_check_intersection2(APoints, BPoints, with_z, 1);
1572 
1573  return ret;
1574 }
#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
GRASS_INTERPFL_EXPORT struct line_pnts * Pnts
GRASS_INTERPFL_EXPORT int count
int Vect_line_check_intersection2(struct line_pnts *APoints, struct line_pnts *BPoints, int with_z)
Check if 2 lines intersect.
Definition: intersect2.c:1542
int line_check_intersection2(struct line_pnts *APoints, struct line_pnts *BPoints, int with_z, int all)
Definition: intersect2.c:1262
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:1564
#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:677
#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