GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
vector/Vlib/intersect.c
Go to the documentation of this file.
1 
21 #include <stdlib.h>
22 #include <math.h>
23 #include <grass/gis.h>
24 #include <grass/Vect.h>
25 #include <grass/glocale.h>
26 
27 /* function prototypes */
28 static int cmp_cross(const void *pa, const void *pb);
29 static void add_cross(int asegment, double adistance, int bsegment,
30  double bdistance, double x, double y);
31 static double dist2(double x1, double y1, double x2, double y2);
32 
33 #if 0
34 static int ident(double x1, double y1, double x2, double y2, double thresh);
35 #endif
36 static int cross_seg(int id, int *arg);
37 static int find_cross(int id, int *arg);
38 
39 
40 /* Some parts of code taken from grass50 v.spag/linecros.c
41  *
42  * Based on the following:
43  *
44  * (ax2-ax1)r1 - (bx2-bx1)r2 = ax2 - ax1
45  * (ay2-ay1)r1 - (by2-by1)r2 = ay2 - ay1
46  *
47  * Solving for r1 and r2, if r1 and r2 are between 0 and 1,
48  * then line segments (ax1,ay1)(ax2,ay2) and (bx1,by1)(bx2,by2)
49  * intersect
50  * ****************************************************************/
51 
52 #define D ((ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2))
53 #define D1 ((bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2))
54 #define D2 ((ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1))
55 
56 /* Intersect 2 line segments.
57  *
58  * Returns: 0 - do not intersect
59  * 1 - intersect at one point
60  * \ / \ / \ /
61  * \/ \/ \/
62  * /\ \
63  * / \ \
64  * 2 - partial overlap ( \/ )
65  * ------ a ( distance < threshold )
66  * ------ b ( )
67  * 3 - a contains b ( /\ )
68  * ---------- a ----------- a
69  * ---- b ----- b
70  * 4 - b contains a
71  * ---- a ----- a
72  * ---------- b ----------- b
73  * 5 - identical
74  * ---------- a
75  * ---------- b
76  *
77  * Intersection points:
78  * return point1 breakes: point2 breaks: distance1 on: distance2 on:
79  * 0 - - - -
80  * 1 a,b - a b
81  * 2 a b a b
82  * 3 a a a a
83  * 4 b b b b
84  * 5 - - - -
85  *
86  * Sometimes (often) is important to get the same coordinates for a x b and b x a.
87  * To reach this, the segments a,b are 'sorted' at the beginning, so that for the same switched segments,
88  * results are identical. (reason is that double values are always rounded because of limited number
89  * of decimal places and for different order of coordinates, the results would be different)
90  *
91  */
92 
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, r2, dtol, t;
117  int switched = 0;
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  /* Check identical lines */
132  if ((ax1 == bx1 && ay1 == by1 && ax2 == bx2 && ay2 == by2) ||
133  (ax1 == bx2 && ay1 == by2 && ax2 == bx1 && ay2 == by1)) {
134  G_debug(2, " -> identical segments");
135  *x1 = ax1;
136  *y1 = ay1;
137  *z1 = az1;
138  *x2 = ax2;
139  *y2 = ay2;
140  *z2 = az2;
141  return 5;
142  }
143 
144  /* 'Sort' lines by x1, x2, y1, y2 */
145  if (bx1 < ax1)
146  switched = 1;
147  else if (bx1 == ax1) {
148  if (bx2 < ax2)
149  switched = 1;
150  else if (bx2 == ax2) {
151  if (by1 < ay1)
152  switched = 1;
153  else if (by1 == ay1) {
154  if (by2 < ay2)
155  switched = 1; /* by2 != ay2 (would be identical */
156  }
157  }
158  }
159  if (switched) {
160  t = ax1;
161  ax1 = bx1;
162  bx1 = t;
163  t = ay1;
164  ay1 = by1;
165  by1 = t;
166  t = ax2;
167  ax2 = bx2;
168  bx2 = t;
169  t = ay2;
170  ay2 = by2;
171  by2 = t;
172  }
173 
174  d = D;
175  d1 = D1;
176  d2 = D2;
177 
178  G_debug(2, "Vect_segment_intersection(): d = %f, d1 = %f, d2 = %f", d, d1,
179  d2);
180 
181  /* TODO: dtol was originaly set to 1.0e-10, which was usualy working but not always.
182  * Can it be a problem to set the tolerance to 0.0 ? */
183  dtol = 0.0;
184  if (fabs(d) > dtol) {
185  r1 = D1 / d;
186  r2 = D2 / d;
187 
188  G_debug(2, " -> not parallel/collinear: r1 = %f, r2 = %f", r1, r2);
189 
190  if (r1 < 0 || r1 > 1 || r2 < 0 || r2 > 1) {
191  G_debug(2, " -> no intersection");
192  return 0;
193  }
194 
195  *x1 = ax1 + r1 * (ax2 - ax1);
196  *y1 = ay1 + r1 * (ay2 - ay1);
197  *z1 = 0;
198 
199  G_debug(2, " -> intersection %f, %f", *x1, *y1);
200  return 1;
201  }
202 
203  /* segments are parallel or collinear */
204  G_debug(3, " -> parallel/collinear");
205 
206  if (D1 || D2) { /* lines are parallel */
207  G_debug(2, " -> parallel");
208  return 0;
209  }
210 
211  /* segments are colinear. check for overlap */
212 
213  /* Collinear vertical */
214  /* original code assumed lines were not both vertical
215  * so there is a special case if they are */
216  if (ax1 == ax2 && bx1 == bx2 && ax1 == bx1) {
217  G_debug(2, " -> collinear vertical");
218  if (ay1 > ay2) {
219  t = ay1;
220  ay1 = ay2;
221  ay2 = t;
222  } /* to be sure that ay1 < ay2 */
223  if (by1 > by2) {
224  t = by1;
225  by1 = by2;
226  by2 = t;
227  } /* to be sure that by1 < by2 */
228  if (ay1 > by2 || ay2 < by1) {
229  G_debug(2, " -> no intersection");
230  return 0;
231  }
232 
233  /* end points */
234  if (ay1 == by2) {
235  *x1 = ax1;
236  *y1 = ay1;
237  *z1 = 0;
238  G_debug(2, " -> connected by end points");
239  return 1; /* endpoints only */
240  }
241  if (ay2 == by1) {
242  *x1 = ax2;
243  *y1 = ay2;
244  *z1 = 0;
245  G_debug(2, " -> connected by end points");
246  return 1; /* endpoints only */
247  }
248 
249  /* heneral overlap */
250  G_debug(3, " -> vertical overlap");
251  /* a contains b */
252  if (ay1 <= by1 && ay2 >= by2) {
253  G_debug(2, " -> a contains b");
254  *x1 = bx1;
255  *y1 = by1;
256  *z1 = 0;
257  *x2 = bx2;
258  *y2 = by2;
259  *z2 = 0;
260  if (!switched)
261  return 3;
262  else
263  return 4;
264  }
265  /* b contains a */
266  if (ay1 >= by1 && ay2 <= by2) {
267  G_debug(2, " -> b contains a");
268  *x1 = ax1;
269  *y1 = ay1;
270  *z1 = 0;
271  *x2 = ax2;
272  *y2 = ay2;
273  *z2 = 0;
274  if (!switched)
275  return 4;
276  else
277  return 3;
278  }
279 
280  /* general overlap, 2 intersection points */
281  G_debug(2, " -> partial overlap");
282  if (by1 > ay1 && by1 < ay2) { /* b1 in a */
283  if (!switched) {
284  *x1 = bx1;
285  *y1 = by1;
286  *z1 = 0;
287  *x2 = ax2;
288  *y2 = ay2;
289  *z2 = 0;
290  }
291  else {
292  *x1 = ax2;
293  *y1 = ay2;
294  *z1 = 0;
295  *x2 = bx1;
296  *y2 = by1;
297  *z2 = 0;
298  }
299  return 2;
300  }
301  if (by2 > ay1 && by2 < ay2) { /* b2 in a */
302  if (!switched) {
303  *x1 = bx2;
304  *y1 = by2;
305  *z1 = 0;
306  *x2 = ax1;
307  *y2 = ay1;
308  *z2 = 0;
309  }
310  else {
311  *x1 = ax1;
312  *y1 = ay1;
313  *z1 = 0;
314  *x2 = bx2;
315  *y2 = by2;
316  *z2 = 0;
317  }
318  return 2;
319  }
320 
321  /* should not be reached */
322  G_warning(_("Vect_segment_intersection() ERROR (collinear vertical segments)"));
323  G_warning("%.15g %.15g", ax1, ay1);
324  G_warning("%.15g %.15g", ax2, ay2);
325  G_warning("x");
326  G_warning("%.15g %.15g", bx1, by1);
327  G_warning("%.15g %.15g", bx2, by2);
328 
329  return 0;
330  }
331 
332  G_debug(2, " -> collinear non vertical");
333 
334  /* Collinear non vertical */
335  if ((bx1 > ax1 && bx2 > ax1 && bx1 > ax2 && bx2 > ax2) ||
336  (bx1 < ax1 && bx2 < ax1 && bx1 < ax2 && bx2 < ax2)) {
337  G_debug(2, " -> no intersection");
338  return 0;
339  }
340 
341  /* there is overlap or connected end points */
342  G_debug(2, " -> overlap/connected end points");
343 
344  /* end points */
345  if ((ax1 == bx1 && ay1 == by1) || (ax1 == bx2 && ay1 == by2)) {
346  *x1 = ax1;
347  *y1 = ay1;
348  *z1 = 0;
349  G_debug(2, " -> connected by end points");
350  return 1;
351  }
352  if ((ax2 == bx1 && ay2 == by1) || (ax2 == bx2 && ay2 == by2)) {
353  *x1 = ax2;
354  *y1 = ay2;
355  *z1 = 0;
356  G_debug(2, " -> connected by end points");
357  return 1;
358  }
359 
360  if (ax1 > ax2) {
361  t = ax1;
362  ax1 = ax2;
363  ax2 = t;
364  t = ay1;
365  ay1 = ay2;
366  ay2 = t;
367  } /* to be sure that ax1 < ax2 */
368  if (bx1 > bx2) {
369  t = bx1;
370  bx1 = bx2;
371  bx2 = t;
372  t = by1;
373  by1 = by2;
374  by2 = t;
375  } /* to be sure that bx1 < bx2 */
376 
377  /* a contains b */
378  if (ax1 <= bx1 && ax2 >= bx2) {
379  G_debug(2, " -> a contains b");
380  *x1 = bx1;
381  *y1 = by1;
382  *z1 = 0;
383  *x2 = bx2;
384  *y2 = by2;
385  *z2 = 0;
386  if (!switched)
387  return 3;
388  else
389  return 4;
390  }
391  /* b contains a */
392  if (ax1 >= bx1 && ax2 <= bx2) {
393  G_debug(2, " -> b contains a");
394  *x1 = ax1;
395  *y1 = ay1;
396  *z1 = 0;
397  *x2 = ax2;
398  *y2 = ay2;
399  *z2 = 0;
400  if (!switched)
401  return 4;
402  else
403  return 3;
404  }
405 
406  /* general overlap, 2 intersection points (lines are not vertical) */
407  G_debug(2, " -> partial overlap");
408  if (bx1 > ax1 && bx1 < ax2) { /* b1 is in a */
409  if (!switched) {
410  *x1 = bx1;
411  *y1 = by1;
412  *z1 = 0;
413  *x2 = ax2;
414  *y2 = ay2;
415  *z2 = 0;
416  }
417  else {
418  *x1 = ax2;
419  *y1 = ay2;
420  *z1 = 0;
421  *x2 = bx1;
422  *y2 = by1;
423  *z2 = 0;
424  }
425  return 2;
426  }
427  if (bx2 > ax1 && bx2 < ax2) { /* b2 is in a */
428  if (!switched) {
429  *x1 = bx2;
430  *y1 = by2;
431  *z1 = 0;
432  *x2 = ax1;
433  *y2 = ay1;
434  *z2 = 0;
435  }
436  else {
437  *x1 = ax1;
438  *y1 = ay1;
439  *z1 = 0;
440  *x2 = bx2;
441  *y2 = by2;
442  *z2 = 0;
443  }
444  return 2;
445  }
446 
447  /* should not be reached */
448  G_warning(_("Vect_segment_intersection() ERROR (collinear non vertical segments)"));
449  G_warning("%.15g %.15g", ax1, ay1);
450  G_warning("%.15g %.15g", ax2, ay2);
451  G_warning("x");
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 whic 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, int *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 = *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 
600 int
601 Vect_line_intersection(struct line_pnts *APoints,
602  struct line_pnts *BPoints,
603  struct line_pnts ***ALines,
604  struct line_pnts ***BLines,
605  int *nalines, int *nblines, int with_z)
606 {
607  int i, j, k, l, last_seg, seg, last;
608  int n_alive_cross;
609  double dist, curdist, last_dist, last_x, last_y, last_z;
610  double x, y, rethresh;
611  struct Rect rect;
612  struct line_pnts **XLines, *Points;
613  struct Node *RTree;
614  struct line_pnts *Points1, *Points2; /* first, second points */
615  int seg1, seg2, vert1, vert2;
616 
617  n_cross = 0;
618  rethresh = 0.000001; /* TODO */
619  APnts = APoints;
620  BPnts = BPoints;
621 
622  /* RE (representation error).
623  * RE thresh above is nonsense of course, the RE threshold should be based on
624  * number of significant digits for double (IEEE-754) which is 15 or 16 and exponent.
625  * The number above is in fact not required threshold, and will not work
626  * for example: equator length is 40.075,695 km (8 digits), units are m (+3)
627  * and we want precision in mm (+ 3) = 14 -> minimum rethresh may be around 0.001
628  * ?Maybe all nonsense? */
629 
630  /* Warning: This function is also used to intersect the line by itself i.e. APoints and
631  * BPoints are identical. I am not sure if it is clever, but it seems to work, but
632  * we have to keep this in mind and handle some special cases (maybe) */
633 
634  /* TODO: 3D, RE threshold, GV_POINTS (line x point) */
635 
636  /* Take each segment from A and intersect by each segment from B.
637  *
638  * All intersections are found first and saved to array, then sorted by a distance along the line,
639  * and then the line is split to pieces.
640  *
641  * Note: If segments are collinear, check if previous/next segments are also collinear,
642  * in that case do not break:
643  * +----------+
644  * +----+-----+ etc.
645  * doesn't need to be broken
646  *
647  * Note: If 2 adjacent segments of line B have common vertex exactly (or within thresh) on line A,
648  * intersection points for these B segments may differ due to RE:
649  * ------------ a ----+--+---- ----+--+----
650  * /\ => / \ or maybe \/
651  * b0 / \ b1 / \ even: /\
652  *
653  * -> solution: snap all breaks to nearest vertices first within RE threshold
654  *
655  * Question: Snap all breaks to each other within RE threshold?
656  *
657  * Note: If a break is snapped to end point or two breaks are snapped to the same vertex
658  * resulting new line is degenerated => before line is added to array, it must be checked
659  * if line is not degenerated
660  *
661  * Note: to snap to vertices is important for cases where line A is broken by B and C line
662  * at the same point:
663  * \ / b no snap \ /
664  * \/ could ----+--+----
665  * ------ a result
666  * /\ in ?: /\
667  * / \ c / \
668  *
669  * Note: once we snap breaks to vertices, we have to do that for both lines A and B in the same way
670  * and because we cannot be sure that A childrens will not change a bit by break(s)
671  * we have to break both A and B at once i.e. in one Vect_line_intersection () call.
672  */
673 
674  /* Spatial index: lines may be very long (thousands of segments) and check each segment
675  * with each from second line takes a long time (n*m). Because of that, spatial index
676  * is build first for the second line and segments from the first line are broken by segments
677  * in bound box */
678 
679  /* Create rtree for B line */
680  RTree = RTreeNewIndex();
681  for (i = 0; i < BPoints->n_points - 1; i++) {
682  if (BPoints->x[i] <= BPoints->x[i + 1]) {
683  rect.boundary[0] = BPoints->x[i];
684  rect.boundary[3] = BPoints->x[i + 1];
685  }
686  else {
687  rect.boundary[0] = BPoints->x[i + 1];
688  rect.boundary[3] = BPoints->x[i];
689  }
690 
691  if (BPoints->y[i] <= BPoints->y[i + 1]) {
692  rect.boundary[1] = BPoints->y[i];
693  rect.boundary[4] = BPoints->y[i + 1];
694  }
695  else {
696  rect.boundary[1] = BPoints->y[i + 1];
697  rect.boundary[4] = BPoints->y[i];
698  }
699 
700  if (BPoints->z[i] <= BPoints->z[i + 1]) {
701  rect.boundary[2] = BPoints->z[i];
702  rect.boundary[5] = BPoints->z[i + 1];
703  }
704  else {
705  rect.boundary[2] = BPoints->z[i + 1];
706  rect.boundary[5] = BPoints->z[i];
707  }
708 
709  RTreeInsertRect(&rect, i + 1, &RTree, 0); /* B line segment numbers in rtree start from 1 */
710  }
711 
712  /* Break segments in A by segments in B */
713  for (i = 0; i < APoints->n_points - 1; i++) {
714  if (APoints->x[i] <= APoints->x[i + 1]) {
715  rect.boundary[0] = APoints->x[i];
716  rect.boundary[3] = APoints->x[i + 1];
717  }
718  else {
719  rect.boundary[0] = APoints->x[i + 1];
720  rect.boundary[3] = APoints->x[i];
721  }
722 
723  if (APoints->y[i] <= APoints->y[i + 1]) {
724  rect.boundary[1] = APoints->y[i];
725  rect.boundary[4] = APoints->y[i + 1];
726  }
727  else {
728  rect.boundary[1] = APoints->y[i + 1];
729  rect.boundary[4] = APoints->y[i];
730  }
731  if (APoints->z[i] <= APoints->z[i + 1]) {
732  rect.boundary[2] = APoints->z[i];
733  rect.boundary[5] = APoints->z[i + 1];
734  }
735  else {
736  rect.boundary[2] = APoints->z[i + 1];
737  rect.boundary[5] = APoints->z[i];
738  }
739 
740  j = RTreeSearch(RTree, &rect, (void *)cross_seg, &i); /* A segment number from 0 */
741  }
742 
743  /* Free RTree */
744  RTreeDestroyNode(RTree);
745 
746  G_debug(2, "n_cross = %d", n_cross);
747  /* Lines do not cross each other */
748  if (n_cross == 0) {
749  *nalines = 0;
750  *nblines = 0;
751  return 0;
752  }
753 
754  /* Snap breaks to nearest vertices within RE threshold */
755  for (i = 0; i < n_cross; i++) {
756  /* 1. of A seg */
757  seg = cross[i].segment[0];
758  curdist =
759  dist2(cross[i].x, cross[i].y, APoints->x[seg], APoints->y[seg]);
760  x = APoints->x[seg];
761  y = APoints->y[seg];
762 
763  /* 2. of A seg */
764  dist =
765  dist2(cross[i].x, cross[i].y, APoints->x[seg + 1],
766  APoints->y[seg + 1]);
767  if (dist < curdist) {
768  curdist = dist;
769  x = APoints->x[seg + 1];
770  y = APoints->y[seg + 1];
771  }
772 
773  /* 1. of B seg */
774  seg = cross[i].segment[1];
775  dist =
776  dist2(cross[i].x, cross[i].y, BPoints->x[seg], BPoints->y[seg]);
777  if (dist < curdist) {
778  curdist = dist;
779  x = BPoints->x[seg];
780  y = BPoints->y[seg];
781  }
782  dist = dist2(cross[i].x, cross[i].y, BPoints->x[seg + 1], BPoints->y[seg + 1]); /* 2. of B seg */
783  if (dist < curdist) {
784  curdist = dist;
785  x = BPoints->x[seg + 1];
786  y = BPoints->y[seg + 1];
787  }
788  if (curdist < rethresh * rethresh) {
789  cross[i].x = x;
790  cross[i].y = y;
791  }
792  }
793 
794  /* Calculate distances along segments */
795  for (i = 0; i < n_cross; i++) {
796  seg = cross[i].segment[0];
797  cross[i].distance[0] =
798  dist2(APoints->x[seg], APoints->y[seg], cross[i].x, cross[i].y);
799  seg = cross[i].segment[1];
800  cross[i].distance[1] =
801  dist2(BPoints->x[seg], BPoints->y[seg], cross[i].x, cross[i].y);
802  }
803 
804  /* l = 1 ~ line A, l = 2 ~ line B */
805  for (l = 1; l < 3; l++) {
806  for (i = 0; i < n_cross; i++)
807  use_cross[i] = 1;
808 
809  /* Create array of lines */
810  XLines = G_malloc((n_cross + 1) * sizeof(struct line_pnts *));
811 
812  if (l == 1) {
813  G_debug(2, "Clean and create array for line A");
814  Points = APoints;
815  Points1 = APoints;
816  Points2 = BPoints;
817  current = 0;
818  second = 1;
819  }
820  else {
821  G_debug(2, "Clean and create array for line B");
822  Points = BPoints;
823  Points1 = BPoints;
824  Points2 = APoints;
825  current = 1;
826  second = 0;
827  }
828 
829  /* Sort points along lines */
830  qsort((void *)cross, sizeof(char) * n_cross, sizeof(CROSS),
831  cmp_cross);
832 
833  /* Print all (raw) breaks */
834  for (i = 0; i < n_cross; i++) {
835  G_debug(3,
836  " cross = %d seg1/dist1 = %d/%f seg2/dist2 = %d/%f x = %f y = %f",
837  i, cross[i].segment[current],
838  sqrt(cross[i].distance[current]),
839  cross[i].segment[second], sqrt(cross[i].distance[second]),
840  cross[i].x, cross[i].y);
841  }
842 
843  /* Remove breaks on first/last line vertices */
844  for (i = 0; i < n_cross; i++) {
845  if (use_cross[i] == 1) {
846  j = Points1->n_points - 1;
847 
848  /* Note: */
849  if ((cross[i].segment[current] == 0 &&
850  cross[i].x == Points1->x[0] &&
851  cross[i].y == Points1->y[0]) ||
852  (cross[i].segment[current] == j - 1 &&
853  cross[i].x == Points1->x[j] &&
854  cross[i].y == Points1->y[j])) {
855  use_cross[i] = 0; /* first/last */
856  G_debug(3, "cross %d deleted (first/last point)", i);
857  }
858  }
859  }
860 
861  /* Remove breaks with collinear previous and next segments on 1 and 2 */
862  /* Note: breaks with collinear previous and nex must be remove duplicates,
863  * otherwise some cross may be lost. Example (+ is vertex):
864  * B first cross intersections: A/B segment:
865  * | 0/0, 0/1, 1/0, 1/1 - collinear previous and next
866  * AB -----+----+--- A 0/4, 0/5, 1/4, 1/5 - OK
867  * \___|
868  * B
869  * This should not inluence that break is always on first segment, see below (I hope)
870  */
871  /* TODO: this doesn't find identical with breaks on revious/next */
872  for (i = 0; i < n_cross; i++) {
873  if (use_cross[i] == 0)
874  continue;
875  G_debug(3, " is %d between colinear?", i);
876 
877  seg1 = cross[i].segment[current];
878  seg2 = cross[i].segment[second];
879 
880  /* Is it vertex on 1, which? */
881  if (cross[i].x == Points1->x[seg1] &&
882  cross[i].y == Points1->y[seg1]) {
883  vert1 = seg1;
884  }
885  else if (cross[i].x == Points1->x[seg1 + 1] &&
886  cross[i].y == Points1->y[seg1 + 1]) {
887  vert1 = seg1 + 1;
888  }
889  else {
890  G_debug(3, " -> is not vertex on 1. line");
891  continue;
892  }
893 
894  /* Is it vertex on 2, which? */
895  /* For 1. line it is easy, because breaks on vertex are always at end vertex
896  * for 2. line we need to find which vertex is on break if any (vert2 starts from 0) */
897  if (cross[i].x == Points2->x[seg2] &&
898  cross[i].y == Points2->y[seg2]) {
899  vert2 = seg2;
900  }
901  else if (cross[i].x == Points2->x[seg2 + 1] &&
902  cross[i].y == Points2->y[seg2 + 1]) {
903  vert2 = seg2 + 1;
904  }
905  else {
906  G_debug(3, " -> is not vertex on 2. line");
907  continue;
908  }
909  G_debug(3, " seg1/vert1 = %d/%d seg2/ver2 = %d/%d", seg1,
910  vert1, seg2, vert2);
911 
912  /* Check if the second vertex is not first/last */
913  if (vert2 == 0 || vert2 == Points2->n_points - 1) {
914  G_debug(3, " -> vertex 2 (%d) is first/last", vert2);
915  continue;
916  }
917 
918  /* Are there first vertices of this segment identical */
919  if (!((Points1->x[vert1 - 1] == Points2->x[vert2 - 1] &&
920  Points1->y[vert1 - 1] == Points2->y[vert2 - 1] &&
921  Points1->x[vert1 + 1] == Points2->x[vert2 + 1] &&
922  Points1->y[vert1 + 1] == Points2->y[vert2 + 1]) ||
923  (Points1->x[vert1 - 1] == Points2->x[vert2 + 1] &&
924  Points1->y[vert1 - 1] == Points2->y[vert2 + 1] &&
925  Points1->x[vert1 + 1] == Points2->x[vert2 - 1] &&
926  Points1->y[vert1 + 1] == Points2->y[vert2 - 1])
927  )
928  ) {
929  G_debug(3, " -> previous/next are not identical");
930  continue;
931  }
932 
933  use_cross[i] = 0;
934 
935  G_debug(3, " -> collinear -> remove");
936  }
937 
938  /* Remove duplicates, i.e. merge all identical breaks to one.
939  * We must be careful because two points with identical coordinates may be distant if measured along
940  * the line:
941  * | Segments b0 and b1 overlap, b0 runs up, b1 down.
942  * | Two inersections may be merged for a, because they are identical,
943  * -----+---- a but cannot be merged for b, because both b0 and b1 must be broken.
944  * | I.e. Breaks on b have identical coordinates, but there are not identical
945  * b0 | b1 if measured along line b.
946  *
947  * -> Breaks may be merged as identical if lay on the same segment, or on vertex connecting
948  * 2 adjacent segments the points lay on
949  *
950  * Note: if duplicate is on a vertex, the break is removed from next segment =>
951  * break on vertex is always on first segment of this vertex (used below)
952  */
953  last = -1;
954  for (i = 1; i < n_cross; i++) {
955  if (use_cross[i] == 0)
956  continue;
957  if (last == -1) { /* set first alive */
958  last = i;
959  continue;
960  }
961  seg = cross[i].segment[current];
962  /* compare with last */
963  G_debug(3, " duplicate ?: cross = %d seg = %d dist = %f", i,
964  cross[i].segment[current], cross[i].distance[current]);
965  if ((cross[i].segment[current] == cross[last].segment[current] &&
966  cross[i].distance[current] == cross[last].distance[current])
967  || (cross[i].segment[current] ==
968  cross[last].segment[current] + 1 &&
969  cross[i].distance[current] == 0 &&
970  cross[i].x == cross[last].x &&
971  cross[i].y == cross[last].y)) {
972  G_debug(3, " cross %d identical to last -> removed", i);
973  use_cross[i] = 0; /* identical */
974  }
975  else {
976  last = i;
977  }
978  }
979 
980  /* Create array of new lines */
981  /* Count alive crosses */
982  n_alive_cross = 0;
983  G_debug(3, " alive crosses:");
984  for (i = 0; i < n_cross; i++) {
985  if (use_cross[i] == 1) {
986  G_debug(3, " %d", i);
987  n_alive_cross++;
988  }
989  }
990 
991  k = 0;
992  if (n_alive_cross > 0) {
993  /* Add last line point at the end of cross array (cross alley) */
994  use_cross[n_cross] = 1;
995  j = Points->n_points - 1;
996  cross[n_cross].x = Points->x[j];
997  cross[n_cross].y = Points->y[j];
998  cross[n_cross].segment[current] = Points->n_points - 2;
999 
1000  last_seg = 0;
1001  last_dist = 0;
1002  last_x = Points->x[0];
1003  last_y = Points->y[0];
1004  last_z = Points->z[0];
1005  /* Go through all cross (+last line point) and create for each new line
1006  * starting at last_* and ending at cross (last point) */
1007  for (i = 0; i <= n_cross; i++) { /* i.e. n_cross + 1 new lines */
1008  seg = cross[i].segment[current];
1009  G_debug(2, "%d seg = %d dist = %f", i, seg,
1010  cross[i].distance[current]);
1011  if (use_cross[i] == 0) {
1012  G_debug(3, " removed -> next");
1013  continue;
1014  }
1015 
1016  G_debug(2, " New line:");
1017  XLines[k] = Vect_new_line_struct();
1018  /* add last intersection or first point first */
1019  Vect_append_point(XLines[k], last_x, last_y, last_z);
1020  G_debug(2, " append last vert: %f %f", last_x, last_y);
1021 
1022  /* add first points of segments between last and current seg */
1023  for (j = last_seg + 1; j <= seg; j++) {
1024  G_debug(2, " segment j = %d", j);
1025  /* skipp vertex identical to last break */
1026  if ((j == last_seg + 1) && Points->x[j] == last_x &&
1027  Points->y[j] == last_y) {
1028  G_debug(2, " -> skip (identical to last break)");
1029  continue;
1030  }
1031  Vect_append_point(XLines[k], Points->x[j], Points->y[j],
1032  Points->z[j]);
1033  G_debug(2, " append first of seg: %f %f", Points->x[j],
1034  Points->y[j]);
1035  }
1036 
1037  /* add current cross or end point */
1038  Vect_append_point(XLines[k], cross[i].x, cross[i].y, 0.0);
1039  G_debug(2, " append cross / last point: %f %f", cross[i].x,
1040  cross[i].y);
1041  last_seg = seg;
1042  last_x = cross[i].x;
1043  last_y = cross[i].y, last_z = 0;
1044 
1045  /* Check if line is degenerate */
1046  if (dig_line_degenerate(XLines[k]) > 0) {
1047  G_debug(2, " line is degenerate -> skipped");
1048  Vect_destroy_line_struct(XLines[k]);
1049  }
1050  else {
1051  k++;
1052  }
1053  }
1054  }
1055  if (l == 1) {
1056  *nalines = k;
1057  *ALines = XLines;
1058  }
1059  else {
1060  *nblines = k;
1061  *BLines = XLines;
1062  }
1063  }
1064 
1065  return 1;
1066 }
1067 
1068 static struct line_pnts *APnts, *BPnts, *IPnts;
1069 
1070 static int cross_found; /* set by find_cross() */
1071 
1072 /* break segments (called by rtree search) */
1073 static int find_cross(int id, int *arg)
1074 {
1075  double x1, y1, z1, x2, y2, z2;
1076  int i, j, ret;
1077 
1078  /* !!! segment number for B lines is returned as +1 */
1079  i = *arg;
1080  j = id - 1;
1081  /* Note: -1 to make up for the +1 when data was inserted */
1082 
1083  ret = Vect_segment_intersection(APnts->x[i], APnts->y[i], APnts->z[i],
1084  APnts->x[i + 1], APnts->y[i + 1],
1085  APnts->z[i + 1], BPnts->x[j], BPnts->y[j],
1086  BPnts->z[j], BPnts->x[j + 1],
1087  BPnts->y[j + 1], BPnts->z[j + 1], &x1,
1088  &y1, &z1, &x2, &y2, &z2, 0);
1089 
1090  if (!IPnts)
1091  IPnts = Vect_new_line_struct();
1092 
1093  switch (ret) {
1094  case 0:
1095  case 5:
1096  break;
1097  case 1:
1098  if (0 > Vect_copy_xyz_to_pnts(IPnts, &x1, &y1, &z1, 1))
1099  G_warning(_("Error while adding point to array. Out of memory"));
1100  break;
1101  case 2:
1102  case 3:
1103  case 4:
1104  if (0 > Vect_copy_xyz_to_pnts(IPnts, &x1, &y1, &z1, 1))
1105  G_warning(_("Error while adding point to array. Out of memory"));
1106  if (0 > Vect_copy_xyz_to_pnts(IPnts, &x2, &y2, &z2, 1))
1107  G_warning(_("Error while adding point to array. Out of memory"));
1108  break;
1109  }
1110  /* add ALL (including end points and duplicates), clean later */
1111  if (ret > 0) {
1112  cross_found = 1;
1113  return 0;
1114  }
1115  return 1; /* keep going */
1116 }
1117 
1130 int
1131 Vect_line_check_intersection(struct line_pnts *APoints,
1132  struct line_pnts *BPoints, int with_z)
1133 {
1134  int i;
1135  double dist, rethresh;
1136  struct Rect rect;
1137  struct Node *RTree;
1138 
1139  rethresh = 0.000001; /* TODO */
1140  APnts = APoints;
1141  BPnts = BPoints;
1142 
1143  /* TODO: 3D, RE (representation error) threshold, GV_POINTS (line x point) */
1144 
1145  if (!IPnts)
1146  IPnts = Vect_new_line_struct();
1147 
1148  /* If one or both are point (Points->n_points == 1) */
1149  if (APoints->n_points == 1 && BPoints->n_points == 1) {
1150  if (APoints->x[0] == BPoints->x[0] && APoints->y[0] == BPoints->y[0]) {
1151  if (!with_z) {
1152  if (0 >
1153  Vect_copy_xyz_to_pnts(IPnts, &APoints->x[0],
1154  &APoints->y[0], NULL, 1))
1155  G_warning(_("Error while adding point to array. Out of memory"));
1156  return 1;
1157  }
1158  else {
1159  if (APoints->z[0] == BPoints->z[0]) {
1160  if (0 >
1161  Vect_copy_xyz_to_pnts(IPnts, &APoints->x[0],
1162  &APoints->y[0], &APoints->z[0],
1163  1))
1164  G_warning(_("Error while adding point to array. Out of memory"));
1165  return 1;
1166  }
1167  else
1168  return 0;
1169  }
1170  }
1171  else {
1172  return 0;
1173  }
1174  }
1175 
1176  if (APoints->n_points == 1) {
1177  Vect_line_distance(BPoints, APoints->x[0], APoints->y[0],
1178  APoints->z[0], with_z, NULL, NULL, NULL, &dist,
1179  NULL, NULL);
1180 
1181  if (dist <= rethresh) {
1182  if (0 >
1183  Vect_copy_xyz_to_pnts(IPnts, &APoints->x[0], &APoints->y[0],
1184  &APoints->z[0], 1))
1185  G_warning(_("Error while adding point to array. Out of memory"));
1186  return 1;
1187  }
1188  else {
1189  return 0;
1190  }
1191  }
1192 
1193  if (BPoints->n_points == 1) {
1194  Vect_line_distance(APoints, BPoints->x[0], BPoints->y[0],
1195  BPoints->z[0], with_z, NULL, NULL, NULL, &dist,
1196  NULL, NULL);
1197 
1198  if (dist <= rethresh) {
1199  if (0 >
1200  Vect_copy_xyz_to_pnts(IPnts, &BPoints->x[0], &BPoints->y[0],
1201  &BPoints->z[0], 1))
1202  G_warning(_("Error while adding point to array. Out of memory"));
1203  return 1;
1204  }
1205  else
1206  return 0;
1207  }
1208 
1209  /* Take each segment from A and find if intersects any segment from B. */
1210 
1211  /* Spatial index: lines may be very long (thousands of segments) and check each segment
1212  * with each from second line takes a long time (n*m). Because of that, spatial index
1213  * is build first for the second line and segments from the first line are broken by segments
1214  * in bound box */
1215 
1216  /* Create rtree for B line */
1217  RTree = RTreeNewIndex();
1218  for (i = 0; i < BPoints->n_points - 1; i++) {
1219  if (BPoints->x[i] <= BPoints->x[i + 1]) {
1220  rect.boundary[0] = BPoints->x[i];
1221  rect.boundary[3] = BPoints->x[i + 1];
1222  }
1223  else {
1224  rect.boundary[0] = BPoints->x[i + 1];
1225  rect.boundary[3] = BPoints->x[i];
1226  }
1227 
1228  if (BPoints->y[i] <= BPoints->y[i + 1]) {
1229  rect.boundary[1] = BPoints->y[i];
1230  rect.boundary[4] = BPoints->y[i + 1];
1231  }
1232  else {
1233  rect.boundary[1] = BPoints->y[i + 1];
1234  rect.boundary[4] = BPoints->y[i];
1235  }
1236 
1237  if (BPoints->z[i] <= BPoints->z[i + 1]) {
1238  rect.boundary[2] = BPoints->z[i];
1239  rect.boundary[5] = BPoints->z[i + 1];
1240  }
1241  else {
1242  rect.boundary[2] = BPoints->z[i + 1];
1243  rect.boundary[5] = BPoints->z[i];
1244  }
1245 
1246  RTreeInsertRect(&rect, i + 1, &RTree, 0); /* B line segment numbers in rtree start from 1 */
1247  }
1248 
1249  /* Find intersection */
1250  cross_found = 0;
1251 
1252  for (i = 0; i < APoints->n_points - 1; i++) {
1253  if (APoints->x[i] <= APoints->x[i + 1]) {
1254  rect.boundary[0] = APoints->x[i];
1255  rect.boundary[3] = APoints->x[i + 1];
1256  }
1257  else {
1258  rect.boundary[0] = APoints->x[i + 1];
1259  rect.boundary[3] = APoints->x[i];
1260  }
1261 
1262  if (APoints->y[i] <= APoints->y[i + 1]) {
1263  rect.boundary[1] = APoints->y[i];
1264  rect.boundary[4] = APoints->y[i + 1];
1265  }
1266  else {
1267  rect.boundary[1] = APoints->y[i + 1];
1268  rect.boundary[4] = APoints->y[i];
1269  }
1270  if (APoints->z[i] <= APoints->z[i + 1]) {
1271  rect.boundary[2] = APoints->z[i];
1272  rect.boundary[5] = APoints->z[i + 1];
1273  }
1274  else {
1275  rect.boundary[2] = APoints->z[i + 1];
1276  rect.boundary[5] = APoints->z[i];
1277  }
1278 
1279  RTreeSearch(RTree, &rect, (void *)find_cross, &i); /* A segment number from 0 */
1280 
1281  if (cross_found) {
1282  break;
1283  }
1284  }
1285 
1286  /* Free RTree */
1287  RTreeDestroyNode(RTree);
1288 
1289  return cross_found;
1290 }
1291 
1305 int
1306 Vect_line_get_intersections(struct line_pnts *APoints,
1307  struct line_pnts *BPoints,
1308  struct line_pnts *IPoints, int with_z)
1309 {
1310  int ret;
1311 
1312  IPnts = IPoints;
1313  ret = Vect_line_check_intersection(APoints, BPoints, with_z);
1314 
1315  return ret;
1316 }