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