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