GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-eff64031f2
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[out] ALines array of new lines created from original A line
652  * \param[out] BLines array of new lines created from original B line
653  * \param[out] nalines number of new lines (ALines)
654  * \param[out] nblines number of new lines (BLines)
655  * \param with_z 3D, not supported!
656  *
657  * \return 0 no intersection
658  * \return 1 intersection found
659  */
660 int Vect_line_intersection(struct line_pnts *APoints, struct line_pnts *BPoints,
661  struct bound_box *ABox, struct bound_box *BBox,
662  struct line_pnts ***ALines,
663  struct line_pnts ***BLines, int *nalines,
664  int *nblines, int with_z UNUSED)
665 {
666  int i, j, k, l, last_seg, seg, last;
667  int n_alive_cross;
668  double dist, curdist, last_x, last_y, last_z;
669  double x, y, rethresh;
670  struct line_pnts **XLines, *Points;
671  struct RTree *MyRTree;
672  struct line_pnts *Points1, *Points2; /* first, second points */
673  int seg1, seg2, vert1, vert2;
674  static struct RTree_Rect rect;
675  static int rect_init = 0;
676  struct bound_box box, abbox;
677 
678  if (debug_level == -1) {
679  const char *dstr = G_getenv_nofatal("DEBUG");
680 
681  if (dstr != NULL)
682  debug_level = atoi(dstr);
683  else
684  debug_level = 0;
685  }
686 
687  if (!rect_init) {
688  rect.boundary = G_malloc(6 * sizeof(RectReal));
689  rect_init = 6;
690  }
691 
692  n_cross = 0;
693  rethresh = 0.000001; /* TODO */
694  APnts = APoints;
695  BPnts = BPoints;
696 
697  /* RE (representation error).
698  * RE thresh above is nonsense of course, the RE threshold should be based
699  * on number of significant digits for double (IEEE-754) which is 15 or 16
700  * and exponent. The number above is in fact not required threshold, and
701  * will not work for example: equator length is 40.075,695 km (8 digits),
702  * units are m (+3) and we want precision in mm (+ 3) = 14 -> minimum
703  * rethresh may be around 0.001 ?Maybe all nonsense? Use rounding error of
704  * the unit in the least place ? max of fabs(x), fabs(y) rethresh = pow(2,
705  * log2(max) - 53) */
706 
707  /* Warning: This function is also used to intersect the line by itself i.e.
708  * APoints and BPoints are identical. I am not sure if it is clever, but it
709  * seems to work, but we have to keep this in mind and handle some special
710  * cases (maybe) */
711 
712  /* TODO: 3D, RE threshold, GV_POINTS (line x point) */
713 
714  /* Take each segment from A and intersect by each segment from B.
715  *
716  * All intersections are found first and saved to array, then sorted by a
717  * distance along the line, and then the line is split to pieces.
718  *
719  * Note: If segments are collinear, check if previous/next segments are
720  * also collinear, in that case do not break:
721  * +----------+
722  * +----+-----+ etc.
723  * doesn't need to be broken
724  *
725  * Note: If 2 adjacent segments of line B have common vertex exactly (or
726  * within thresh) on line A, intersection points for these B segments may
727  * differ due to RE:
728  * ------------ a ----+--+---- ----+--+----
729  * /\ => / \ or maybe \/
730  * b0 / \ b1 / \ even: /\
731  *
732  * -> solution: snap all breaks to nearest vertices first within RE
733  * threshold
734  *
735  * Question: Snap all breaks to each other within RE threshold?
736  *
737  * Note: If a break is snapped to end point or two breaks are snapped to
738  * the same vertex resulting new line is degenerated => before line is added
739  * to array, it must be checked if line is not degenerated
740  *
741  * Note: to snap to vertices is important for cases where line A is broken
742  * by B and C line at the same point: \ / b no snap \ /
743  * \/ could ----+--+----
744  * ------ a result
745  * /\ in ?: /\
746  * / \ c / \
747  *
748  * Note: once we snap breaks to vertices, we have to do that for both lines
749  * A and B in the same way and because we cannot be sure that A children
750  * will not change a bit by break(s) we have to break both A and B at once
751  * i.e. in one Vect_line_intersection () call.
752  */
753 
754  /* Spatial index: lines may be very long (thousands of segments) and check
755  * each segment with each from second line takes a long time (n*m). Because
756  * of that, spatial index is build first for the second line and segments
757  * from the first line are broken by segments in bound box */
758 
759  if (!Vect_box_overlap(ABox, BBox)) {
760  *nalines = 0;
761  *nblines = 0;
762  return 0;
763  }
764 
765  abbox = *BBox;
766  if (abbox.N > ABox->N)
767  abbox.N = ABox->N;
768  if (abbox.S < ABox->S)
769  abbox.S = ABox->S;
770  if (abbox.E > ABox->E)
771  abbox.E = ABox->E;
772  if (abbox.W < ABox->W)
773  abbox.W = ABox->W;
774  if (abbox.T > ABox->T)
775  abbox.T = ABox->T;
776  if (abbox.B < ABox->B)
777  abbox.B = ABox->B;
778 
779  abbox.N += rethresh;
780  abbox.S -= rethresh;
781  abbox.E += rethresh;
782  abbox.W -= rethresh;
783  abbox.T += rethresh;
784  abbox.B -= rethresh;
785 
786  /* Create rtree for B line */
787  MyRTree = RTreeCreateTree(-1, 0, 2);
788  RTreeSetOverflow(MyRTree, 0);
789  for (i = 0; i < BPoints->n_points - 1; i++) {
790  if (BPoints->x[i] <= BPoints->x[i + 1]) {
791  rect.boundary[0] = BPoints->x[i];
792  rect.boundary[3] = BPoints->x[i + 1];
793  }
794  else {
795  rect.boundary[0] = BPoints->x[i + 1];
796  rect.boundary[3] = BPoints->x[i];
797  }
798 
799  if (BPoints->y[i] <= BPoints->y[i + 1]) {
800  rect.boundary[1] = BPoints->y[i];
801  rect.boundary[4] = BPoints->y[i + 1];
802  }
803  else {
804  rect.boundary[1] = BPoints->y[i + 1];
805  rect.boundary[4] = BPoints->y[i];
806  }
807 
808  if (BPoints->z[i] <= BPoints->z[i + 1]) {
809  rect.boundary[2] = BPoints->z[i];
810  rect.boundary[5] = BPoints->z[i + 1];
811  }
812  else {
813  rect.boundary[2] = BPoints->z[i + 1];
814  rect.boundary[5] = BPoints->z[i];
815  }
816 
817  box.W = rect.boundary[0] - rethresh;
818  box.S = rect.boundary[1] - rethresh;
819  box.B = rect.boundary[2] - rethresh;
820  box.E = rect.boundary[3] + rethresh;
821  box.N = rect.boundary[4] + rethresh;
822  box.T = rect.boundary[5] + rethresh;
823 
824  if (Vect_box_overlap(&abbox, &box)) {
826  &rect, i + 1,
827  MyRTree); /* B line segment numbers in rtree start from 1 */
828  }
829  }
830 
831  /* Break segments in A by segments in B */
832  for (i = 0; i < APoints->n_points - 1; i++) {
833  if (APoints->x[i] <= APoints->x[i + 1]) {
834  rect.boundary[0] = APoints->x[i];
835  rect.boundary[3] = APoints->x[i + 1];
836  }
837  else {
838  rect.boundary[0] = APoints->x[i + 1];
839  rect.boundary[3] = APoints->x[i];
840  }
841 
842  if (APoints->y[i] <= APoints->y[i + 1]) {
843  rect.boundary[1] = APoints->y[i];
844  rect.boundary[4] = APoints->y[i + 1];
845  }
846  else {
847  rect.boundary[1] = APoints->y[i + 1];
848  rect.boundary[4] = APoints->y[i];
849  }
850  if (APoints->z[i] <= APoints->z[i + 1]) {
851  rect.boundary[2] = APoints->z[i];
852  rect.boundary[5] = APoints->z[i + 1];
853  }
854  else {
855  rect.boundary[2] = APoints->z[i + 1];
856  rect.boundary[5] = APoints->z[i];
857  }
858  box.W = rect.boundary[0] - rethresh;
859  box.S = rect.boundary[1] - rethresh;
860  box.B = rect.boundary[2] - rethresh;
861  box.E = rect.boundary[3] + rethresh;
862  box.N = rect.boundary[4] + rethresh;
863  box.T = rect.boundary[5] + rethresh;
864 
865  if (Vect_box_overlap(&abbox, &box)) {
866  j = RTreeSearch(MyRTree, &rect, cross_seg,
867  &i); /* A segment number from 0 */
868  }
869  }
870 
871  /* Free RTree */
872  RTreeDestroyTree(MyRTree);
873 
874  G_debug(2, "n_cross = %d", n_cross);
875  /* Lines do not cross each other */
876  if (n_cross == 0) {
877  *nalines = 0;
878  *nblines = 0;
879  return 0;
880  }
881 
882  /* Snap breaks to nearest vertices within RE threshold */
883  /* Calculate distances along segments */
884  for (i = 0; i < n_cross; i++) {
885 
886  /* 1. of A seg */
887  seg = cross[i].segment[0];
888  curdist =
889  dist2(cross[i].x, cross[i].y, APoints->x[seg], APoints->y[seg]);
890  x = APoints->x[seg];
891  y = APoints->y[seg];
892 
893  cross[i].distance[0] = curdist;
894 
895  /* 2. of A seg */
896  dist = dist2(cross[i].x, cross[i].y, APoints->x[seg + 1],
897  APoints->y[seg + 1]);
898  if (dist < curdist) {
899  curdist = dist;
900  x = APoints->x[seg + 1];
901  y = APoints->y[seg + 1];
902  }
903 
904  /* 1. of B seg */
905  seg = cross[i].segment[1];
906  dist = dist2(cross[i].x, cross[i].y, BPoints->x[seg], BPoints->y[seg]);
907  cross[i].distance[1] = dist;
908 
909  if (dist < curdist) {
910  curdist = dist;
911  x = BPoints->x[seg];
912  y = BPoints->y[seg];
913  }
914  /* 2. of B seg */
915  dist = dist2(cross[i].x, cross[i].y, BPoints->x[seg + 1],
916  BPoints->y[seg + 1]);
917  if (dist < curdist) {
918  curdist = dist;
919  x = BPoints->x[seg + 1];
920  y = BPoints->y[seg + 1];
921  }
922  if (curdist < rethresh * rethresh) {
923  cross[i].x = x;
924  cross[i].y = y;
925 
926  /* Update distances along segments */
927  seg = cross[i].segment[0];
928  cross[i].distance[0] =
929  dist2(APoints->x[seg], APoints->y[seg], cross[i].x, cross[i].y);
930  seg = cross[i].segment[1];
931  cross[i].distance[1] =
932  dist2(BPoints->x[seg], BPoints->y[seg], cross[i].x, cross[i].y);
933  }
934  }
935 
936  /* l = 1 ~ line A, l = 2 ~ line B */
937  for (l = 1; l < 3; l++) {
938  for (i = 0; i < n_cross; i++)
939  use_cross[i] = 1;
940 
941  /* Create array of lines */
942  XLines = G_malloc((n_cross + 1) * sizeof(struct line_pnts *));
943 
944  if (l == 1) {
945  G_debug(2, "Clean and create array for line A");
946  Points = APoints;
947  Points1 = APoints;
948  Points2 = BPoints;
949  current = 0;
950  second = 1;
951  }
952  else {
953  G_debug(2, "Clean and create array for line B");
954  Points = BPoints;
955  Points1 = BPoints;
956  Points2 = APoints;
957  current = 1;
958  second = 0;
959  }
960 
961  /* Sort points along lines */
962  qsort((void *)cross, sizeof(char) * n_cross, sizeof(CROSS), cmp_cross);
963 
964  /* Print all (raw) breaks */
965  /* avoid loop when not debugging */
966  if (debug_level > 2) {
967  for (i = 0; i < n_cross; i++) {
968  G_debug(
969  3,
970  " cross = %d seg1/dist1 = %d/%f seg2/dist2 = %d/%f x = %f "
971  "y = %f",
972  i, cross[i].segment[current],
973  sqrt(cross[i].distance[current]), cross[i].segment[second],
974  sqrt(cross[i].distance[second]), cross[i].x, cross[i].y);
975  }
976  }
977 
978  /* Remove breaks on first/last line vertices */
979  for (i = 0; i < n_cross; i++) {
980  if (use_cross[i] == 1) {
981  j = Points1->n_points - 1;
982 
983  /* Note: */
984  if ((cross[i].segment[current] == 0 &&
985  cross[i].x == Points1->x[0] &&
986  cross[i].y == Points1->y[0]) ||
987  (cross[i].segment[current] == j - 1 &&
988  cross[i].x == Points1->x[j] &&
989  cross[i].y == Points1->y[j])) {
990  use_cross[i] = 0; /* first/last */
991  G_debug(3, "cross %d deleted (first/last point)", i);
992  }
993  }
994  }
995 
996  /* Remove breaks with collinear previous and next segments on 1 and 2 */
997  /* Note: breaks with collinear previous and nex must be remove
998  * duplicates, otherwise some cross may be lost. Example (+ is vertex):
999  * B first cross intersections: A/B segment:
1000  * | 0/0, 0/1, 1/0, 1/1 - collinear previous
1001  * and next AB -----+----+--- A 0/4, 0/5, 1/4, 1/5 - OK
1002  * \___|
1003  * B
1004  * This should not influence that break is always on first segment, see
1005  * below (I hope)
1006  */
1007  /* TODO: this doesn't find identical with breaks on revious/next */
1008  for (i = 0; i < n_cross; i++) {
1009  if (use_cross[i] == 0)
1010  continue;
1011  G_debug(3, " is %d between colinear?", i);
1012 
1013  seg1 = cross[i].segment[current];
1014  seg2 = cross[i].segment[second];
1015 
1016  /* Is it vertex on 1, which? */
1017  if (cross[i].x == Points1->x[seg1] &&
1018  cross[i].y == Points1->y[seg1]) {
1019  vert1 = seg1;
1020  }
1021  else if (cross[i].x == Points1->x[seg1 + 1] &&
1022  cross[i].y == Points1->y[seg1 + 1]) {
1023  vert1 = seg1 + 1;
1024  }
1025  else {
1026  G_debug(3, " -> is not vertex on 1. line");
1027  continue;
1028  }
1029 
1030  /* Is it vertex on 2, which? */
1031  /* For 1. line it is easy, because breaks on vertex are always at
1032  * end vertex for 2. line we need to find which vertex is on break
1033  * if any (vert2 starts from 0) */
1034  if (cross[i].x == Points2->x[seg2] &&
1035  cross[i].y == Points2->y[seg2]) {
1036  vert2 = seg2;
1037  }
1038  else if (cross[i].x == Points2->x[seg2 + 1] &&
1039  cross[i].y == Points2->y[seg2 + 1]) {
1040  vert2 = seg2 + 1;
1041  }
1042  else {
1043  G_debug(3, " -> is not vertex on 2. line");
1044  continue;
1045  }
1046  G_debug(3, " seg1/vert1 = %d/%d seg2/ver2 = %d/%d", seg1, vert1,
1047  seg2, vert2);
1048 
1049  /* Check if the second vertex is not first/last */
1050  if (vert2 == 0 || vert2 == Points2->n_points - 1) {
1051  G_debug(3, " -> vertex 2 (%d) is first/last", vert2);
1052  continue;
1053  }
1054 
1055  /* Are there first vertices of this segment identical */
1056  if (!((Points1->x[vert1 - 1] == Points2->x[vert2 - 1] &&
1057  Points1->y[vert1 - 1] == Points2->y[vert2 - 1] &&
1058  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  G_debug(3, " -> previous/next are not identical");
1065  continue;
1066  }
1067 
1068  use_cross[i] = 0;
1069 
1070  G_debug(3, " -> collinear -> remove");
1071  }
1072 
1073  /* Remove duplicates, i.e. merge all identical breaks to one.
1074  * We must be careful because two points with identical coordinates may
1075  * be distant if measured along the line: | Segments b0 and b1
1076  * overlap, b0 runs up, b1 down. | Two inersections may be
1077  * merged for a, because they are identical,
1078  * -----+---- a but cannot be merged for b, because both b0 and b1
1079  * must be broken. | I.e. Breaks on b have identical
1080  * coordinates, but there are not identical b0 | b1 if measured
1081  * along line b.
1082  *
1083  * -> Breaks may be merged as identical if lay on the same segment,
1084  * or on vertex connecting 2 adjacent segments the points lay on
1085  *
1086  * Note: if duplicate is on a vertex, the break is removed from next
1087  * segment => break on vertex is always on first segment of this vertex
1088  * (used below)
1089  */
1090  last = -1;
1091  for (i = 0; i < n_cross; i++) {
1092  if (use_cross[i] == 0)
1093  continue;
1094  if (last == -1) { /* set first alive */
1095  last = i;
1096  continue;
1097  }
1098  seg = cross[i].segment[current];
1099  /* compare with last */
1100  G_debug(3, " duplicate ?: cross = %d seg = %d dist = %f", i,
1101  cross[i].segment[current], cross[i].distance[current]);
1102  if ((cross[i].segment[current] == cross[last].segment[current] &&
1103  cross[i].distance[current] == cross[last].distance[current]) ||
1104  (cross[i].segment[current] ==
1105  cross[last].segment[current] + 1 &&
1106  cross[i].distance[current] == 0 &&
1107  cross[i].x == cross[last].x && cross[i].y == cross[last].y)) {
1108  G_debug(3, " cross %d identical to last -> removed", i);
1109  use_cross[i] = 0; /* identical */
1110  }
1111  else {
1112  last = i;
1113  }
1114  }
1115 
1116  /* Create array of new lines */
1117  /* Count alive crosses */
1118  n_alive_cross = 0;
1119  G_debug(3, " alive crosses:");
1120  for (i = 0; i < n_cross; i++) {
1121  if (use_cross[i] == 1) {
1122  G_debug(3, " %d", i);
1123  n_alive_cross++;
1124  }
1125  }
1126 
1127  k = 0;
1128  if (n_alive_cross > 0) {
1129  /* Add last line point at the end of cross array (cross alley) */
1130  use_cross[n_cross] = 1;
1131  j = Points->n_points - 1;
1132  cross[n_cross].x = Points->x[j];
1133  cross[n_cross].y = Points->y[j];
1134  cross[n_cross].segment[current] = Points->n_points - 2;
1135 
1136  last_seg = 0;
1137  last_x = Points->x[0];
1138  last_y = Points->y[0];
1139  last_z = Points->z[0];
1140  /* Go through all cross (+last line point) and create for each new
1141  * line starting at last_* and ending at cross (last point) */
1142  for (i = 0; i <= n_cross; i++) { /* i.e. n_cross + 1 new lines */
1143  seg = cross[i].segment[current];
1144  G_debug(2, "%d seg = %d dist = %f", i, seg,
1145  cross[i].distance[current]);
1146  if (use_cross[i] == 0) {
1147  G_debug(3, " removed -> next");
1148  continue;
1149  }
1150 
1151  G_debug(2, " New line:");
1152  XLines[k] = Vect_new_line_struct();
1153  /* add last intersection or first point first */
1154  Vect_append_point(XLines[k], last_x, last_y, last_z);
1155  G_debug(2, " append last vert: %f %f", last_x, last_y);
1156 
1157  /* add first points of segments between last and current seg */
1158  for (j = last_seg + 1; j <= seg; j++) {
1159  G_debug(2, " segment j = %d", j);
1160  /* skip vertex identical to last break */
1161  if ((j == last_seg + 1) && Points->x[j] == last_x &&
1162  Points->y[j] == last_y) {
1163  G_debug(2, " -> skip (identical to last break)");
1164  continue;
1165  }
1166  Vect_append_point(XLines[k], Points->x[j], Points->y[j],
1167  Points->z[j]);
1168  G_debug(2, " append first of seg: %f %f", Points->x[j],
1169  Points->y[j]);
1170  }
1171 
1172  last_seg = seg;
1173  last_x = cross[i].x;
1174  last_y = cross[i].y;
1175  last_z = 0;
1176  /* calculate last_z */
1177  if (Points->z[last_seg] == Points->z[last_seg + 1]) {
1178  last_z = Points->z[last_seg + 1];
1179  }
1180  else if (last_x == Points->x[last_seg] &&
1181  last_y == Points->y[last_seg]) {
1182  last_z = Points->z[last_seg];
1183  }
1184  else if (last_x == Points->x[last_seg + 1] &&
1185  last_y == Points->y[last_seg + 1]) {
1186  last_z = Points->z[last_seg + 1];
1187  }
1188  else {
1189  dist = dist2(Points->x[last_seg], Points->x[last_seg + 1],
1190  Points->y[last_seg], Points->y[last_seg + 1]);
1191  last_z =
1192  (Points->z[last_seg] *
1193  sqrt(cross[i].distance[current]) +
1194  Points->z[last_seg + 1] *
1195  (sqrt(dist) - sqrt(cross[i].distance[current]))) /
1196  sqrt(dist);
1197  }
1198 
1199  /* add current cross or end point */
1200  Vect_append_point(XLines[k], cross[i].x, cross[i].y, last_z);
1201  G_debug(2, " append cross / last point: %f %f", cross[i].x,
1202  cross[i].y);
1203 
1204  /* Check if line is degenerate */
1205  if (dig_line_degenerate(XLines[k]) > 0) {
1206  G_debug(2, " line is degenerate -> skipped");
1207  Vect_destroy_line_struct(XLines[k]);
1208  }
1209  else {
1210  k++;
1211  }
1212  }
1213  }
1214  if (l == 1) {
1215  *nalines = k;
1216  *ALines = XLines;
1217  }
1218  else {
1219  *nblines = k;
1220  *BLines = XLines;
1221  }
1222  }
1223 
1224  /* clean up */
1225 
1226  return 1;
1227 }
1228 
1229 static struct line_pnts *APnts, *BPnts, *IPnts;
1230 
1231 static int cross_found; /* set by find_cross() */
1232 static int report_all; /* should all crossings be reported or just first one */
1233 
1234 /* break segments (called by rtree search) */
1235 static int find_cross(int id, const struct RTree_Rect *rect UNUSED, void *arg)
1236 {
1237  double x1, y1, z1, x2, y2, z2;
1238  int i, j, ret;
1239 
1240  /* !!! segment number for B lines is returned as +1 */
1241  i = *(int *)arg;
1242  j = id - 1;
1243  /* Note: -1 to make up for the +1 when data was inserted */
1244 
1246  APnts->x[i], APnts->y[i], APnts->z[i], APnts->x[i + 1], APnts->y[i + 1],
1247  APnts->z[i + 1], BPnts->x[j], BPnts->y[j], BPnts->z[j], BPnts->x[j + 1],
1248  BPnts->y[j + 1], BPnts->z[j + 1], &x1, &y1, &z1, &x2, &y2, &z2, 0);
1249 
1250  switch (ret) {
1251  case 0:
1252  case 5:
1253  break;
1254  case 1:
1255  if (0 > Vect_append_point(IPnts, x1, y1, z1))
1256  G_warning(_("Error while adding point to array. Out of memory"));
1257  break;
1258  case 2:
1259  case 3:
1260  case 4:
1261  if (0 > Vect_append_point(IPnts, x1, y1, z1))
1262  G_warning(_("Error while adding point to array. Out of memory"));
1263  if (0 > Vect_append_point(IPnts, x2, y2, z2))
1264  G_warning(_("Error while adding point to array. Out of memory"));
1265  break;
1266  }
1267  /* add ALL (including end points and duplicates), clean later */
1268  if (ret > 0) {
1269  cross_found = 1;
1270  if (!report_all)
1271  return 0;
1272  }
1273  return 1; /* keep going */
1274 }
1275 
1277  struct line_pnts *BPoints, int with_z)
1278 {
1279  int i;
1280  double dist, rethresh;
1281  struct RTree *MyRTree;
1282  static struct RTree_Rect rect;
1283  static int rect_init = 0;
1284 
1285  if (!rect_init) {
1286  rect.boundary = G_malloc(6 * sizeof(RectReal));
1287  rect_init = 6;
1288  }
1289 
1290  rethresh = 0.000001; /* TODO */
1291  APnts = APoints;
1292  BPnts = BPoints;
1293 
1294  /* TODO: 3D, RE (representation error) threshold, GV_POINTS (line x point)
1295  */
1296 
1297  if (!IPnts)
1298  IPnts = Vect_new_line_struct();
1299  Vect_reset_line(IPnts);
1300 
1301  /* If one or both are point (Points->n_points == 1) */
1302  if (APoints->n_points == 1 && BPoints->n_points == 1) {
1303  if (APoints->x[0] == BPoints->x[0] && APoints->y[0] == BPoints->y[0]) {
1304  if (!with_z) {
1305  if (report_all &&
1306  0 > Vect_copy_xyz_to_pnts(IPnts, &APoints->x[0],
1307  &APoints->y[0], NULL, 1))
1308  G_warning(
1309  _("Error while adding point to array. Out of memory"));
1310  return 1;
1311  }
1312  else {
1313  if (APoints->z[0] == BPoints->z[0]) {
1314  if (report_all &&
1315  0 > Vect_copy_xyz_to_pnts(IPnts, &APoints->x[0],
1316  &APoints->y[0],
1317  &APoints->z[0], 1))
1318  G_warning(_("Error while adding point to array. Out of "
1319  "memory"));
1320  return 1;
1321  }
1322  else
1323  return 0;
1324  }
1325  }
1326  else {
1327  return 0;
1328  }
1329  }
1330 
1331  if (APoints->n_points == 1) {
1332  Vect_line_distance(BPoints, APoints->x[0], APoints->y[0], APoints->z[0],
1333  with_z, NULL, NULL, NULL, &dist, NULL, NULL);
1334 
1335  if (dist <= rethresh) {
1336  if (report_all &&
1337  0 > Vect_copy_xyz_to_pnts(IPnts, &APoints->x[0], &APoints->y[0],
1338  &APoints->z[0], 1))
1339  G_warning(
1340  _("Error while adding point to array. Out of memory"));
1341  return 1;
1342  }
1343  else {
1344  return 0;
1345  }
1346  }
1347 
1348  if (BPoints->n_points == 1) {
1349  Vect_line_distance(APoints, BPoints->x[0], BPoints->y[0], BPoints->z[0],
1350  with_z, NULL, NULL, NULL, &dist, NULL, NULL);
1351 
1352  if (dist <= rethresh) {
1353  if (report_all &&
1354  0 > Vect_copy_xyz_to_pnts(IPnts, &BPoints->x[0], &BPoints->y[0],
1355  &BPoints->z[0], 1))
1356  G_warning(
1357  _("Error while adding point to array. Out of memory"));
1358  return 1;
1359  }
1360  else
1361  return 0;
1362  }
1363 
1364  /* Take each segment from A and find if intersects any segment from B. */
1365 
1366  /* Spatial index: lines may be very long (thousands of segments) and check
1367  * each segment with each from second line takes a long time (n*m). Because
1368  * of that, spatial index is build first for the second line and segments
1369  * from the first line are broken by segments in bound box */
1370 
1371  /* Create rtree for B line */
1372  MyRTree = RTreeCreateTree(-1, 0, 2);
1373  RTreeSetOverflow(MyRTree, 0);
1374  for (i = 0; i < BPoints->n_points - 1; i++) {
1375  if (BPoints->x[i] <= BPoints->x[i + 1]) {
1376  rect.boundary[0] = BPoints->x[i];
1377  rect.boundary[3] = BPoints->x[i + 1];
1378  }
1379  else {
1380  rect.boundary[0] = BPoints->x[i + 1];
1381  rect.boundary[3] = BPoints->x[i];
1382  }
1383 
1384  if (BPoints->y[i] <= BPoints->y[i + 1]) {
1385  rect.boundary[1] = BPoints->y[i];
1386  rect.boundary[4] = BPoints->y[i + 1];
1387  }
1388  else {
1389  rect.boundary[1] = BPoints->y[i + 1];
1390  rect.boundary[4] = BPoints->y[i];
1391  }
1392 
1393  if (BPoints->z[i] <= BPoints->z[i + 1]) {
1394  rect.boundary[2] = BPoints->z[i];
1395  rect.boundary[5] = BPoints->z[i + 1];
1396  }
1397  else {
1398  rect.boundary[2] = BPoints->z[i + 1];
1399  rect.boundary[5] = BPoints->z[i];
1400  }
1401 
1403  &rect, i + 1,
1404  MyRTree); /* B line segment numbers in rtree start from 1 */
1405  }
1406 
1407  /* Find intersection */
1408  cross_found = 0;
1409 
1410  for (i = 0; i < APoints->n_points - 1; i++) {
1411  if (APoints->x[i] <= APoints->x[i + 1]) {
1412  rect.boundary[0] = APoints->x[i];
1413  rect.boundary[3] = APoints->x[i + 1];
1414  }
1415  else {
1416  rect.boundary[0] = APoints->x[i + 1];
1417  rect.boundary[3] = APoints->x[i];
1418  }
1419 
1420  if (APoints->y[i] <= APoints->y[i + 1]) {
1421  rect.boundary[1] = APoints->y[i];
1422  rect.boundary[4] = APoints->y[i + 1];
1423  }
1424  else {
1425  rect.boundary[1] = APoints->y[i + 1];
1426  rect.boundary[4] = APoints->y[i];
1427  }
1428  if (APoints->z[i] <= APoints->z[i + 1]) {
1429  rect.boundary[2] = APoints->z[i];
1430  rect.boundary[5] = APoints->z[i + 1];
1431  }
1432  else {
1433  rect.boundary[2] = APoints->z[i + 1];
1434  rect.boundary[5] = APoints->z[i];
1435  }
1436 
1437  RTreeSearch(MyRTree, &rect, find_cross,
1438  &i); /* A segment number from 0 */
1439 
1440  if (!report_all && cross_found) {
1441  break;
1442  }
1443  }
1444 
1445  /* Free RTree */
1446  RTreeDestroyTree(MyRTree);
1447 
1448  return cross_found;
1449 }
1450 
1451 /*!
1452  * \brief Check if 2 lines intersect.
1453  *
1454  * Points (Points->n_points == 1) are also supported.
1455  *
1456  * Superseded by the faster Vect_line_check_intersection2()
1457  * Kept as reference implementation
1458  *
1459  * \param APoints first input line
1460  * \param BPoints second input line
1461  * \param with_z 3D, not supported (only if one or both are points)!
1462  *
1463  * \return 0 no intersection
1464  * \return 1 intersection found
1465  */
1467  struct line_pnts *BPoints, int with_z)
1468 {
1469  report_all = 0;
1470  return line_check_intersection(APoints, BPoints, with_z);
1471 }
1472 
1473 /*!
1474  * \brief Get 2 lines intersection points.
1475  *
1476  * A wrapper around Vect_line_check_intersection() function.
1477  *
1478  * Superseded by the faster Vect_line_get_intersections2()
1479  * Kept as reference implementation
1480  *
1481  * \param APoints first input line
1482  * \param BPoints second input line
1483  * \param[out] IPoints output with intersection points
1484  * \param with_z 3D, not supported (only if one or both are points)!
1485  *
1486  * \return 0 no intersection
1487  * \return 1 intersection found
1488  */
1490  struct line_pnts *BPoints,
1491  struct line_pnts *IPoints, int with_z)
1492 {
1493  int ret;
1494 
1495  report_all = 1;
1496  IPnts = IPoints;
1497  ret = line_check_intersection(APoints, BPoints, with_z);
1498 
1499  return ret;
1500 }
#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:47
#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