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