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