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