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