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