GRASS Programmer's Manual  6.5.svn(2014)-r66266
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);
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, 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  /* Check identical segments */
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 x, y
145  * MUST happen before D, D1, D2 are calculated */
146  switched = 0;
147  if (bx2 < bx1)
148  switched = 1;
149  else if (bx2 == bx1) {
150  if (by2 < by1)
151  switched = 1;
152  }
153
154  if (switched) {
155  t = bx1;
156  bx1 = bx2;
157  bx2 = t;
158  t = by1;
159  by1 = by2;
160  by2 = t;
161  t = bz1;
162  bz1 = bz2;
163  bz2 = t;
164  }
165
166  switched = 0;
167  if (ax2 < ax1)
168  switched = 1;
169  else if (ax2 == ax1) {
170  if (ay2 < ay1)
171  switched = 1;
172  }
173
174  if (switched) {
175  t = ax1;
176  ax1 = ax2;
177  ax2 = t;
178  t = ay1;
179  ay1 = ay2;
180  ay2 = t;
181  t = az1;
182  az1 = az2;
183  az2 = t;
184  }
185
186  d = D;
187  d1 = D1;
188  d2 = D2;
189
190  G_debug(2, "Vect_segment_intersection(): d = %f, d1 = %f, d2 = %f", d, d1,
191  d2);
192
193  /* TODO: dtol was originaly set to 1.0e-10, which was usualy working but not always.
194  * Can it be a problem to set the tolerance to 0.0 ? */
195  dtol = 0.0;
196  if (fabs(d) > dtol) {
197
198  G_debug(2, " -> not parallel/collinear: d1 = %f, d2 = %f", d1, d2);
199  if (d > 0) {
200  if (d1 < 0 || d1 > d || d2 < 0 || d2 > d) {
201  G_debug(2, " -> no intersection");
202  return 0;
203  }
204  }
205  else {
206  if (d1 < d || d1 > 0 || d2 < d || d2 > 0) {
207  G_debug(2, " -> no intersection");
208  return 0;
209  }
210  }
211
212  r1 = d1 / d;
213
214  *x1 = ax1 + r1 * (ax2 - ax1);
215  *y1 = ay1 + r1 * (ay2 - ay1);
216  *z1 = 0;
217
218  G_debug(2, " -> intersection %f, %f", *x1, *y1);
219  return 1;
220  }
221
222  /* segments are parallel or collinear */
223  G_debug(3, " -> parallel/collinear");
224
225  if (d1 || d2) { /* lines are parallel */
226  G_debug(2, " -> parallel");
227  return 0;
228  }
229
230  /* segments are colinear. check for overlap */
231
232  /* Collinear vertical */
233  /* original code assumed lines were not both vertical
234  * so there is a special case if they are */
235  if (ax1 == ax2) {
236  G_debug(2, " -> collinear vertical");
237  if (ay1 > by2 || ay2 < by1) {
238  G_debug(2, " -> no intersection");
239  return 0;
240  }
241
242  /* end points */
243  if (ay1 == by2) {
244  G_debug(2, " -> connected by end points");
245  *x1 = ax1;
246  *y1 = ay1;
247  *z1 = 0;
248
249  return 1; /* endpoints only */
250  }
251  if (ay2 == by1) {
252  G_debug(2, " -> connected by end points");
253  *x1 = ax2;
254  *y1 = ay2;
255  *z1 = 0;
256
257  return 1; /* endpoints only */
258  }
259
260  /* general overlap */
261  G_debug(3, " -> vertical overlap");
262  /* a contains b */
263  if (ay1 <= by1 && ay2 >= by2) {
264  G_debug(2, " -> a contains b");
265  *x1 = bx1;
266  *y1 = by1;
267  *z1 = 0;
268  *x2 = bx2;
269  *y2 = by2;
270  *z2 = 0;
271
272  return 3;
273  }
274  /* b contains a */
275  if (ay1 >= by1 && ay2 <= by2) {
276  G_debug(2, " -> b contains a");
277  *x1 = ax1;
278  *y1 = ay1;
279  *z1 = 0;
280  *x2 = ax2;
281  *y2 = ay2;
282  *z2 = 0;
283
284  return 4;
285  }
286
287  /* general overlap, 2 intersection points */
288  G_debug(2, " -> partial overlap");
289  if (by1 > ay1 && by1 < ay2) { /* b1 in a */
290  G_debug(2, " -> b1 in a");
291  *x1 = bx1;
292  *y1 = by1;
293  *z1 = 0;
294  *x2 = ax2;
295  *y2 = ay2;
296  *z2 = 0;
297
298  return 2;
299  }
300  if (by2 > ay1 && by2 < ay2) { /* b2 in a */
301  G_debug(2, " -> b2 in a");
302  *x1 = ax1;
303  *y1 = ay1;
304  *z1 = 0;
305  *x2 = bx2;
306  *y2 = by2;
307  *z2 = 0;
308
309  return 2;
310  }
311
312  /* should not be reached */
313  G_warning(_("Vect_segment_intersection() ERROR (collinear vertical segments)"));
314  G_warning("a");
315  G_warning("%.15g %.15g", ax1, ay1);
316  G_warning("%.15g %.15g", ax2, ay2);
317  G_warning("b");
318  G_warning("%.15g %.15g", bx1, by1);
319  G_warning("%.15g %.15g", bx2, by2);
320
321  return 0;
322  }
323
324  /* Collinear non vertical */
325
326  G_debug(2, " -> collinear non vertical");
327
328  /* b is to the left or right of a */
329  if ((bx1 > ax2) || (bx2 < ax1)) {
330  /* should not happen if segments are selected from rtree */
331  G_debug(2, " -> no intersection");
332  return 0;
333  }
334
335  /* there is overlap or connected end points */
336  G_debug(2, " -> overlap/connected end points");
337
338  /* end points */
339  if (ax1 == bx2 && ay1 == by2) {
340  G_debug(2, " -> connected by end points");
341  *x1 = ax1;
342  *y1 = ay1;
343  *z1 = 0;
344
345  return 1;
346  }
347  if (ax2 == bx1 && ay2 == by1) {
348  G_debug(2, " -> connected by end points");
349  *x1 = ax2;
350  *y1 = ay2;
351  *z1 = 0;
352
353  return 1;
354  }
355
356  /* a contains b */
357  if (ax1 <= bx1 && ax2 >= bx2) {
358  G_debug(2, " -> a contains b");
359  *x1 = bx1;
360  *y1 = by1;
361  *z1 = 0;
362  *x2 = bx2;
363  *y2 = by2;
364  *z2 = 0;
365
366  return 3;
367  }
368  /* b contains a */
369  if (ax1 >= bx1 && ax2 <= bx2) {
370  G_debug(2, " -> b contains a");
371  *x1 = ax1;
372  *y1 = ay1;
373  *z1 = 0;
374  *x2 = ax2;
375  *y2 = ay2;
376  *z2 = 0;
377
378  return 4;
379  }
380
381  /* general overlap, 2 intersection points (lines are not vertical) */
382  G_debug(2, " -> partial overlap");
383  if (bx1 > ax1 && bx1 < ax2) { /* b1 is in a */
384  G_debug(2, " -> b1 in a");
385  *x1 = bx1;
386  *y1 = by1;
387  *z1 = 0;
388  *x2 = ax2;
389  *y2 = ay2;
390  *z2 = 0;
391
392  return 2;
393  }
394  if (bx2 > ax1 && bx2 < ax2) { /* b2 is in a */
395  G_debug(2, " -> b2 in a");
396  *x1 = ax1;
397  *y1 = ay1;
398  *z1 = 0;
399  *x2 = bx2;
400  *y2 = by2;
401  *z2 = 0;
402
403  return 2;
404  }
405
406  /* should not be reached */
407  G_warning(_("Vect_segment_intersection() ERROR (collinear non vertical segments)"));
408  G_warning("a");
409  G_warning("%.15g %.15g", ax1, ay1);
410  G_warning("%.15g %.15g", ax2, ay2);
411  G_warning("b");
412  G_warning("%.15g %.15g", bx1, by1);
413  G_warning("%.15g %.15g", bx2, by2);
414
415  return 0;
416 }
417
418 typedef struct
419 { /* in arrays 0 - A line , 1 - B line */
420  int segment; /* segment number, start from 0 for first */
421  double distance;
422  double x, y, z;
423 } CROSS;
424
425 /* Current line in arrays is for some functions like cmp() set by: */
426 static int current;
427 static int second; /* line whic is not current */
428
429 static int a_cross = 0;
430 static int n_cross;
431 static CROSS *cross = NULL;
432 static int *use_cross = NULL;
433
435  double bdistance, double x, double y)
436 {
437  if (n_cross == a_cross) {
438  /* Must be space + 1, used later for last line point, do it better */
439  cross =
440  (CROSS *) G_realloc((void *)cross,
441  (a_cross + 101) * sizeof(CROSS));
442  use_cross =
443  (int *)G_realloc((void *)use_cross,
444  (a_cross + 101) * sizeof(int));
445  a_cross += 100;
446  }
447
448  G_debug(5,
449  " add new cross: aseg/dist = %d/%f bseg/dist = %d/%f, x = %f y = %f",
450  asegment, adistance, bsegment, bdistance, x, y);
451  cross[n_cross].segment = asegment;
453  cross[n_cross].segment = bsegment;
454  cross[n_cross].distance = bdistance;
455  cross[n_cross].x = x;
456  cross[n_cross].y = y;
457  n_cross++;
458 }
459
460 static int cmp_cross(const void *pa, const void *pb)
461 {
462  CROSS *p1 = (CROSS *) pa;
463  CROSS *p2 = (CROSS *) pb;
464
465  if (p1->segment[current] < p2->segment[current])
466  return -1;
467  if (p1->segment[current] > p2->segment[current])
468  return 1;
469  /* the same segment */
470  if (p1->distance[current] < p2->distance[current])
471  return -1;
472  if (p1->distance[current] > p2->distance[current])
473  return 1;
474  return 0;
475 }
476
477 static double dist2(double x1, double y1, double x2, double y2)
478 {
479  double dx, dy;
480
481  dx = x2 - x1;
482  dy = y2 - y1;
483  return (dx * dx + dy * dy);
484 }
485
486 #if 0
487 /* returns 1 if points are identical */
488 static int ident(double x1, double y1, double x2, double y2, double thresh)
489 {
490  double dx, dy;
491
492  dx = x2 - x1;
493  dy = y2 - y1;
494  if ((dx * dx + dy * dy) <= thresh * thresh)
495  return 1;
496
497  return 0;
498 }
499 #endif
500
501 /* shared by Vect_line_intersection, Vect_line_check_intersection, cross_seg, find_cross */
502 static struct line_pnts *APnts, *BPnts;
503
504 /* break segments (called by rtree search) */
505 static int cross_seg(int id, int *arg)
506 {
507  double x1, y1, z1, x2, y2, z2;
508  int i, j, ret;
509
510  /* !!! segment number for B lines is returned as +1 */
511  i = *arg;
512  j = id - 1;
513  /* Note: -1 to make up for the +1 when data was inserted */
514
515  ret = Vect_segment_intersection(APnts->x[i], APnts->y[i], APnts->z[i],
516  APnts->x[i + 1], APnts->y[i + 1],
517  APnts->z[i + 1], BPnts->x[j], BPnts->y[j],
518  BPnts->z[j], BPnts->x[j + 1],
519  BPnts->y[j + 1], BPnts->z[j + 1], &x1,
520  &y1, &z1, &x2, &y2, &z2, 0);
521
522  /* add ALL (including end points and duplicates), clean later */
523  if (ret > 0) {
524  G_debug(2, " -> %d x %d: intersection type = %d", i, j, ret);
525  if (ret == 1) { /* one intersection on segment A */
526  G_debug(3, " in %f, %f ", x1, y1);
527  add_cross(i, 0.0, j, 0.0, x1, y1);
528  }
529  else if (ret == 2 || ret == 3 || ret == 4 || ret == 5) {
530  /* partial overlap; a broken in one, b broken in one
531  * or a contains b; a is broken in 2 points (but 1 may be end)
532  * or b contains a; b is broken in 2 points (but 1 may be end)
533  * or identical */
534  G_debug(3, " in %f, %f; %f, %f", x1, y1, x2, y2);
535  add_cross(i, 0.0, j, 0.0, x1, y1);
536  add_cross(i, 0.0, j, 0.0, x2, y2);
537  }
538  }
539  return 1; /* keep going */
540 }
541
560 int
561 Vect_line_intersection(struct line_pnts *APoints,
562  struct line_pnts *BPoints,
563  struct line_pnts ***ALines,
564  struct line_pnts ***BLines,
565  int *nalines, int *nblines, int with_z)
566 {
567  int i, j, k, l, last_seg, seg, last;
568  int n_alive_cross;
569  double dist, curdist, last_dist, last_x, last_y, last_z;
570  double x, y, rethresh;
571  struct Rect rect;
572  struct line_pnts **XLines, *Points;
573  struct Node *RTree;
574  struct line_pnts *Points1, *Points2; /* first, second points */
575  int seg1, seg2, vert1, vert2;
576
577  n_cross = 0;
578  rethresh = 0.000001; /* TODO */
579  APnts = APoints;
580  BPnts = BPoints;
581
582  /* RE (representation error).
583  * RE thresh above is nonsense of course, the RE threshold should be based on
584  * number of significant digits for double (IEEE-754) which is 15 or 16 and exponent.
585  * The number above is in fact not required threshold, and will not work
586  * for example: equator length is 40.075,695 km (8 digits), units are m (+3)
587  * and we want precision in mm (+ 3) = 14 -> minimum rethresh may be around 0.001
588  * ?Maybe all nonsense? */
589
590  /* Warning: This function is also used to intersect the line by itself i.e. APoints and
591  * BPoints are identical. I am not sure if it is clever, but it seems to work, but
592  * we have to keep this in mind and handle some special cases (maybe) */
593
594  /* TODO: 3D, RE threshold, GV_POINTS (line x point) */
595
596  /* Take each segment from A and intersect by each segment from B.
597  *
598  * All intersections are found first and saved to array, then sorted by a distance along the line,
599  * and then the line is split to pieces.
600  *
601  * Note: If segments are collinear, check if previous/next segments are also collinear,
602  * in that case do not break:
603  * +----------+
604  * +----+-----+ etc.
605  * doesn't need to be broken
606  *
607  * Note: If 2 adjacent segments of line B have common vertex exactly (or within thresh) on line A,
608  * intersection points for these B segments may differ due to RE:
609  * ------------ a ----+--+---- ----+--+----
610  * /\ => / \ or maybe \/
611  * b0 / \ b1 / \ even: /\
612  *
613  * -> solution: snap all breaks to nearest vertices first within RE threshold
614  *
615  * Question: Snap all breaks to each other within RE threshold?
616  *
617  * Note: If a break is snapped to end point or two breaks are snapped to the same vertex
618  * resulting new line is degenerated => before line is added to array, it must be checked
619  * if line is not degenerated
620  *
621  * Note: to snap to vertices is important for cases where line A is broken by B and C line
622  * at the same point:
623  * \ / b no snap \ /
624  * \/ could ----+--+----
625  * ------ a result
626  * /\ in ?: /\
627  * / \ c / \
628  *
629  * Note: once we snap breaks to vertices, we have to do that for both lines A and B in the same way
630  * and because we cannot be sure that A childrens will not change a bit by break(s)
631  * we have to break both A and B at once i.e. in one Vect_line_intersection () call.
632  */
633
634  /* Spatial index: lines may be very long (thousands of segments) and check each segment
635  * with each from second line takes a long time (n*m). Because of that, spatial index
636  * is build first for the second line and segments from the first line are broken by segments
637  * in bound box */
638
639  /* Create rtree for B line */
640  RTree = RTreeNewIndex();
641  for (i = 0; i < BPoints->n_points - 1; i++) {
642  if (BPoints->x[i] <= BPoints->x[i + 1]) {
643  rect.boundary = BPoints->x[i];
644  rect.boundary = BPoints->x[i + 1];
645  }
646  else {
647  rect.boundary = BPoints->x[i + 1];
648  rect.boundary = BPoints->x[i];
649  }
650
651  if (BPoints->y[i] <= BPoints->y[i + 1]) {
652  rect.boundary = BPoints->y[i];
653  rect.boundary = BPoints->y[i + 1];
654  }
655  else {
656  rect.boundary = BPoints->y[i + 1];
657  rect.boundary = BPoints->y[i];
658  }
659
660  if (BPoints->z[i] <= BPoints->z[i + 1]) {
661  rect.boundary = BPoints->z[i];
662  rect.boundary = BPoints->z[i + 1];
663  }
664  else {
665  rect.boundary = BPoints->z[i + 1];
666  rect.boundary = BPoints->z[i];
667  }
668
669  RTreeInsertRect(&rect, i + 1, &RTree, 0); /* B line segment numbers in rtree start from 1 */
670  }
671
672  /* Break segments in A by segments in B */
673  for (i = 0; i < APoints->n_points - 1; i++) {
674  if (APoints->x[i] <= APoints->x[i + 1]) {
675  rect.boundary = APoints->x[i];
676  rect.boundary = APoints->x[i + 1];
677  }
678  else {
679  rect.boundary = APoints->x[i + 1];
680  rect.boundary = APoints->x[i];
681  }
682
683  if (APoints->y[i] <= APoints->y[i + 1]) {
684  rect.boundary = APoints->y[i];
685  rect.boundary = APoints->y[i + 1];
686  }
687  else {
688  rect.boundary = APoints->y[i + 1];
689  rect.boundary = APoints->y[i];
690  }
691  if (APoints->z[i] <= APoints->z[i + 1]) {
692  rect.boundary = APoints->z[i];
693  rect.boundary = APoints->z[i + 1];
694  }
695  else {
696  rect.boundary = APoints->z[i + 1];
697  rect.boundary = APoints->z[i];
698  }
699
700  j = RTreeSearch(RTree, &rect, (void *)cross_seg, &i); /* A segment number from 0 */
701  }
702
703  /* Free RTree */
704  RTreeDestroyNode(RTree);
705
706  G_debug(2, "n_cross = %d", n_cross);
707  /* Lines do not cross each other */
708  if (n_cross == 0) {
709  *nalines = 0;
710  *nblines = 0;
711  return 0;
712  }
713
714  /* Snap breaks to nearest vertices within RE threshold */
715  for (i = 0; i < n_cross; i++) {
716  /* 1. of A seg */
717  seg = cross[i].segment;
718  curdist =
719  dist2(cross[i].x, cross[i].y, APoints->x[seg], APoints->y[seg]);
720  x = APoints->x[seg];
721  y = APoints->y[seg];
722
723  /* 2. of A seg */
724  dist =
725  dist2(cross[i].x, cross[i].y, APoints->x[seg + 1],
726  APoints->y[seg + 1]);
727  if (dist < curdist) {
728  curdist = dist;
729  x = APoints->x[seg + 1];
730  y = APoints->y[seg + 1];
731  }
732
733  /* 1. of B seg */
734  seg = cross[i].segment;
735  dist =
736  dist2(cross[i].x, cross[i].y, BPoints->x[seg], BPoints->y[seg]);
737  if (dist < curdist) {
738  curdist = dist;
739  x = BPoints->x[seg];
740  y = BPoints->y[seg];
741  }
742  dist = dist2(cross[i].x, cross[i].y, BPoints->x[seg + 1], BPoints->y[seg + 1]); /* 2. of B seg */
743  if (dist < curdist) {
744  curdist = dist;
745  x = BPoints->x[seg + 1];
746  y = BPoints->y[seg + 1];
747  }
748  if (curdist < rethresh * rethresh) {
749  cross[i].x = x;
750  cross[i].y = y;
751  }
752  }
753
754  /* Calculate distances along segments */
755  for (i = 0; i < n_cross; i++) {
756  seg = cross[i].segment;
757  cross[i].distance =
758  dist2(APoints->x[seg], APoints->y[seg], cross[i].x, cross[i].y);
759  seg = cross[i].segment;
760  cross[i].distance =
761  dist2(BPoints->x[seg], BPoints->y[seg], cross[i].x, cross[i].y);
762  }
763
764  /* l = 1 ~ line A, l = 2 ~ line B */
765  for (l = 1; l < 3; l++) {
766  for (i = 0; i < n_cross; i++)
767  use_cross[i] = 1;
768
769  /* Create array of lines */
770  XLines = G_malloc((n_cross + 1) * sizeof(struct line_pnts *));
771
772  if (l == 1) {
773  G_debug(2, "Clean and create array for line A");
774  Points = APoints;
775  Points1 = APoints;
776  Points2 = BPoints;
777  current = 0;
778  second = 1;
779  }
780  else {
781  G_debug(2, "Clean and create array for line B");
782  Points = BPoints;
783  Points1 = BPoints;
784  Points2 = APoints;
785  current = 1;
786  second = 0;
787  }
788
789  /* Sort points along lines */
790  qsort((void *)cross, sizeof(char) * n_cross, sizeof(CROSS),
791  cmp_cross);
792
793  /* Print all (raw) breaks */
794  for (i = 0; i < n_cross; i++) {
795  G_debug(3,
796  " cross = %d seg1/dist1 = %d/%f seg2/dist2 = %d/%f x = %f y = %f",
797  i, cross[i].segment[current],
798  sqrt(cross[i].distance[current]),
799  cross[i].segment[second], sqrt(cross[i].distance[second]),
800  cross[i].x, cross[i].y);
801  }
802
803  /* Remove breaks on first/last line vertices */
804  for (i = 0; i < n_cross; i++) {
805  if (use_cross[i] == 1) {
806  j = Points1->n_points - 1;
807
808  /* Note: */
809  if ((cross[i].segment[current] == 0 &&
810  cross[i].x == Points1->x &&
811  cross[i].y == Points1->y) ||
812  (cross[i].segment[current] == j - 1 &&
813  cross[i].x == Points1->x[j] &&
814  cross[i].y == Points1->y[j])) {
815  use_cross[i] = 0; /* first/last */
816  G_debug(3, "cross %d deleted (first/last point)", i);
817  }
818  }
819  }
820
821  /* Remove breaks with collinear previous and next segments on 1 and 2 */
822  /* Note: breaks with collinear previous and nex must be remove duplicates,
823  * otherwise some cross may be lost. Example (+ is vertex):
824  * B first cross intersections: A/B segment:
825  * | 0/0, 0/1, 1/0, 1/1 - collinear previous and next
826  * AB -----+----+--- A 0/4, 0/5, 1/4, 1/5 - OK
827  * \___|
828  * B
829  * This should not inluence that break is always on first segment, see below (I hope)
830  */
831  /* TODO: this doesn't find identical with breaks on revious/next */
832  for (i = 0; i < n_cross; i++) {
833  if (use_cross[i] == 0)
834  continue;
835  G_debug(3, " is %d between colinear?", i);
836
837  seg1 = cross[i].segment[current];
838  seg2 = cross[i].segment[second];
839
840  /* Is it vertex on 1, which? */
841  if (cross[i].x == Points1->x[seg1] &&
842  cross[i].y == Points1->y[seg1]) {
843  vert1 = seg1;
844  }
845  else if (cross[i].x == Points1->x[seg1 + 1] &&
846  cross[i].y == Points1->y[seg1 + 1]) {
847  vert1 = seg1 + 1;
848  }
849  else {
850  G_debug(3, " -> is not vertex on 1. line");
851  continue;
852  }
853
854  /* Is it vertex on 2, which? */
855  /* For 1. line it is easy, because breaks on vertex are always at end vertex
856  * for 2. line we need to find which vertex is on break if any (vert2 starts from 0) */
857  if (cross[i].x == Points2->x[seg2] &&
858  cross[i].y == Points2->y[seg2]) {
859  vert2 = seg2;
860  }
861  else if (cross[i].x == Points2->x[seg2 + 1] &&
862  cross[i].y == Points2->y[seg2 + 1]) {
863  vert2 = seg2 + 1;
864  }
865  else {
866  G_debug(3, " -> is not vertex on 2. line");
867  continue;
868  }
869  G_debug(3, " seg1/vert1 = %d/%d seg2/ver2 = %d/%d", seg1,
870  vert1, seg2, vert2);
871
872  /* Check if the second vertex is not first/last */
873  if (vert2 == 0 || vert2 == Points2->n_points - 1) {
874  G_debug(3, " -> vertex 2 (%d) is first/last", vert2);
875  continue;
876  }
877
878  /* Are there first vertices of this segment identical */
879  if (!((Points1->x[vert1 - 1] == Points2->x[vert2 - 1] &&
880  Points1->y[vert1 - 1] == Points2->y[vert2 - 1] &&
881  Points1->x[vert1 + 1] == Points2->x[vert2 + 1] &&
882  Points1->y[vert1 + 1] == Points2->y[vert2 + 1]) ||
883  (Points1->x[vert1 - 1] == Points2->x[vert2 + 1] &&
884  Points1->y[vert1 - 1] == Points2->y[vert2 + 1] &&
885  Points1->x[vert1 + 1] == Points2->x[vert2 - 1] &&
886  Points1->y[vert1 + 1] == Points2->y[vert2 - 1])
887  )
888  ) {
889  G_debug(3, " -> previous/next are not identical");
890  continue;
891  }
892
893  use_cross[i] = 0;
894
895  G_debug(3, " -> collinear -> remove");
896  }
897
898  /* Remove duplicates, i.e. merge all identical breaks to one.
899  * We must be careful because two points with identical coordinates may be distant if measured along
900  * the line:
901  * | Segments b0 and b1 overlap, b0 runs up, b1 down.
902  * | Two inersections may be merged for a, because they are identical,
903  * -----+---- a but cannot be merged for b, because both b0 and b1 must be broken.
904  * | I.e. Breaks on b have identical coordinates, but there are not identical
905  * b0 | b1 if measured along line b.
906  *
907  * -> Breaks may be merged as identical if lay on the same segment, or on vertex connecting
908  * 2 adjacent segments the points lay on
909  *
910  * Note: if duplicate is on a vertex, the break is removed from next segment =>
911  * break on vertex is always on first segment of this vertex (used below)
912  */
913  last = -1;
914  for (i = 1; i < n_cross; i++) {
915  if (use_cross[i] == 0)
916  continue;
917  if (last == -1) { /* set first alive */
918  last = i;
919  continue;
920  }
921  seg = cross[i].segment[current];
922  /* compare with last */
923  G_debug(3, " duplicate ?: cross = %d seg = %d dist = %f", i,
924  cross[i].segment[current], cross[i].distance[current]);
925  if ((cross[i].segment[current] == cross[last].segment[current] &&
926  cross[i].distance[current] == cross[last].distance[current])
927  || (cross[i].segment[current] ==
928  cross[last].segment[current] + 1 &&
929  cross[i].distance[current] == 0 &&
930  cross[i].x == cross[last].x &&
931  cross[i].y == cross[last].y)) {
932  G_debug(3, " cross %d identical to last -> removed", i);
933  use_cross[i] = 0; /* identical */
934  }
935  else {
936  last = i;
937  }
938  }
939
940  /* Create array of new lines */
941  /* Count alive crosses */
942  n_alive_cross = 0;
943  G_debug(3, " alive crosses:");
944  for (i = 0; i < n_cross; i++) {
945  if (use_cross[i] == 1) {
946  G_debug(3, " %d", i);
947  n_alive_cross++;
948  }
949  }
950
951  k = 0;
952  if (n_alive_cross > 0) {
953  /* Add last line point at the end of cross array (cross alley) */
954  use_cross[n_cross] = 1;
955  j = Points->n_points - 1;
956  cross[n_cross].x = Points->x[j];
957  cross[n_cross].y = Points->y[j];
958  cross[n_cross].segment[current] = Points->n_points - 2;
959
960  last_seg = 0;
961  last_dist = 0;
962  last_x = Points->x;
963  last_y = Points->y;
964  last_z = Points->z;
965  /* Go through all cross (+last line point) and create for each new line
966  * starting at last_* and ending at cross (last point) */
967  for (i = 0; i <= n_cross; i++) { /* i.e. n_cross + 1 new lines */
968  seg = cross[i].segment[current];
969  G_debug(2, "%d seg = %d dist = %f", i, seg,
970  cross[i].distance[current]);
971  if (use_cross[i] == 0) {
972  G_debug(3, " removed -> next");
973  continue;
974  }
975
976  G_debug(2, " New line:");
977  XLines[k] = Vect_new_line_struct();
978  /* add last intersection or first point first */
979  Vect_append_point(XLines[k], last_x, last_y, last_z);
980  G_debug(2, " append last vert: %f %f", last_x, last_y);
981
982  /* add first points of segments between last and current seg */
983  for (j = last_seg + 1; j <= seg; j++) {
984  G_debug(2, " segment j = %d", j);
985  /* skipp vertex identical to last break */
986  if ((j == last_seg + 1) && Points->x[j] == last_x &&
987  Points->y[j] == last_y) {
988  G_debug(2, " -> skip (identical to last break)");
989  continue;
990  }
991  Vect_append_point(XLines[k], Points->x[j], Points->y[j],
992  Points->z[j]);
993  G_debug(2, " append first of seg: %f %f", Points->x[j],
994  Points->y[j]);
995  }
996
997  /* add current cross or end point */
998  Vect_append_point(XLines[k], cross[i].x, cross[i].y, 0.0);
999  G_debug(2, " append cross / last point: %f %f", cross[i].x,
1000  cross[i].y);
1001  last_seg = seg;
1002  last_x = cross[i].x;
1003  last_y = cross[i].y, last_z = 0;
1004
1005  /* Check if line is degenerate */
1006  if (dig_line_degenerate(XLines[k]) > 0) {
1007  G_debug(2, " line is degenerate -> skipped");
1008  Vect_destroy_line_struct(XLines[k]);
1009  }
1010  else {
1011  k++;
1012  }
1013  }
1014  }
1015  if (l == 1) {
1016  *nalines = k;
1017  *ALines = XLines;
1018  }
1019  else {
1020  *nblines = k;
1021  *BLines = XLines;
1022  }
1023  }
1024
1025  return 1;
1026 }
1027
1028 static struct line_pnts *APnts, *BPnts, *IPnts;
1029
1030 static int cross_found; /* set by find_cross() */
1031
1032 /* break segments (called by rtree search) */
1033 static int find_cross(int id, int *arg)
1034 {
1035  double x1, y1, z1, x2, y2, z2;
1036  int i, j, ret;
1037
1038  /* !!! segment number for B lines is returned as +1 */
1039  i = *arg;
1040  j = id - 1;
1041  /* Note: -1 to make up for the +1 when data was inserted */
1042
1043  ret = Vect_segment_intersection(APnts->x[i], APnts->y[i], APnts->z[i],
1044  APnts->x[i + 1], APnts->y[i + 1],
1045  APnts->z[i + 1], BPnts->x[j], BPnts->y[j],
1046  BPnts->z[j], BPnts->x[j + 1],
1047  BPnts->y[j + 1], BPnts->z[j + 1], &x1,
1048  &y1, &z1, &x2, &y2, &z2, 0);
1049
1050  if (!IPnts)
1051  IPnts = Vect_new_line_struct();
1052
1053  switch (ret) {
1054  case 0:
1055  case 5:
1056  break;
1057  case 1:
1058  if (0 > Vect_copy_xyz_to_pnts(IPnts, &x1, &y1, &z1, 1))
1059  G_warning(_("Error while adding point to array. Out of memory"));
1060  break;
1061  case 2:
1062  case 3:
1063  case 4:
1064  if (0 > Vect_copy_xyz_to_pnts(IPnts, &x1, &y1, &z1, 1))
1065  G_warning(_("Error while adding point to array. Out of memory"));
1066  if (0 > Vect_copy_xyz_to_pnts(IPnts, &x2, &y2, &z2, 1))
1067  G_warning(_("Error while adding point to array. Out of memory"));
1068  break;
1069  }
1070  /* add ALL (including end points and duplicates), clean later */
1071  if (ret > 0) {
1072  cross_found = 1;
1073  return 0;
1074  }
1075  return 1; /* keep going */
1076 }
1077
1090 int
1091 Vect_line_check_intersection(struct line_pnts *APoints,
1092  struct line_pnts *BPoints, int with_z)
1093 {
1094  int i;
1095  double dist, rethresh;
1096  struct Rect rect;
1097  struct Node *RTree;
1098
1099  rethresh = 0.000001; /* TODO */
1100  APnts = APoints;
1101  BPnts = BPoints;
1102
1103  /* TODO: 3D, RE (representation error) threshold, GV_POINTS (line x point) */
1104
1105  if (!IPnts)
1106  IPnts = Vect_new_line_struct();
1107
1108  /* If one or both are point (Points->n_points == 1) */
1109  if (APoints->n_points == 1 && BPoints->n_points == 1) {
1110  if (APoints->x == BPoints->x && APoints->y == BPoints->y) {
1111  if (!with_z) {
1112  if (0 >
1113  Vect_copy_xyz_to_pnts(IPnts, &APoints->x,
1114  &APoints->y, NULL, 1))
1115  G_warning(_("Error while adding point to array. Out of memory"));
1116  return 1;
1117  }
1118  else {
1119  if (APoints->z == BPoints->z) {
1120  if (0 >
1121  Vect_copy_xyz_to_pnts(IPnts, &APoints->x,
1122  &APoints->y, &APoints->z,
1123  1))
1124  G_warning(_("Error while adding point to array. Out of memory"));
1125  return 1;
1126  }
1127  else
1128  return 0;
1129  }
1130  }
1131  else {
1132  return 0;
1133  }
1134  }
1135
1136  if (APoints->n_points == 1) {
1137  Vect_line_distance(BPoints, APoints->x, APoints->y,
1138  APoints->z, with_z, NULL, NULL, NULL, &dist,
1139  NULL, NULL);
1140
1141  if (dist <= rethresh) {
1142  if (0 >
1143  Vect_copy_xyz_to_pnts(IPnts, &APoints->x, &APoints->y,
1144  &APoints->z, 1))
1145  G_warning(_("Error while adding point to array. Out of memory"));
1146  return 1;
1147  }
1148  else {
1149  return 0;
1150  }
1151  }
1152
1153  if (BPoints->n_points == 1) {
1154  Vect_line_distance(APoints, BPoints->x, BPoints->y,
1155  BPoints->z, with_z, NULL, NULL, NULL, &dist,
1156  NULL, NULL);
1157
1158  if (dist <= rethresh) {
1159  if (0 >
1160  Vect_copy_xyz_to_pnts(IPnts, &BPoints->x, &BPoints->y,
1161  &BPoints->z, 1))
1162  G_warning(_("Error while adding point to array. Out of memory"));
1163  return 1;
1164  }
1165  else
1166  return 0;
1167  }
1168
1169  /* Take each segment from A and find if intersects any segment from B. */
1170
1171  /* Spatial index: lines may be very long (thousands of segments) and check each segment
1172  * with each from second line takes a long time (n*m). Because of that, spatial index
1173  * is build first for the second line and segments from the first line are broken by segments
1174  * in bound box */
1175
1176  /* Create rtree for B line */
1177  RTree = RTreeNewIndex();
1178  for (i = 0; i < BPoints->n_points - 1; i++) {
1179  if (BPoints->x[i] <= BPoints->x[i + 1]) {
1180  rect.boundary = BPoints->x[i];
1181  rect.boundary = BPoints->x[i + 1];
1182  }
1183  else {
1184  rect.boundary = BPoints->x[i + 1];
1185  rect.boundary = BPoints->x[i];
1186  }
1187
1188  if (BPoints->y[i] <= BPoints->y[i + 1]) {
1189  rect.boundary = BPoints->y[i];
1190  rect.boundary = BPoints->y[i + 1];
1191  }
1192  else {
1193  rect.boundary = BPoints->y[i + 1];
1194  rect.boundary = BPoints->y[i];
1195  }
1196
1197  if (BPoints->z[i] <= BPoints->z[i + 1]) {
1198  rect.boundary = BPoints->z[i];
1199  rect.boundary = BPoints->z[i + 1];
1200  }
1201  else {
1202  rect.boundary = BPoints->z[i + 1];
1203  rect.boundary = BPoints->z[i];
1204  }
1205
1206  RTreeInsertRect(&rect, i + 1, &RTree, 0); /* B line segment numbers in rtree start from 1 */
1207  }
1208
1209  /* Find intersection */
1210  cross_found = 0;
1211
1212  for (i = 0; i < APoints->n_points - 1; i++) {
1213  if (APoints->x[i] <= APoints->x[i + 1]) {
1214  rect.boundary = APoints->x[i];
1215  rect.boundary = APoints->x[i + 1];
1216  }
1217  else {
1218  rect.boundary = APoints->x[i + 1];
1219  rect.boundary = APoints->x[i];
1220  }
1221
1222  if (APoints->y[i] <= APoints->y[i + 1]) {
1223  rect.boundary = APoints->y[i];
1224  rect.boundary = APoints->y[i + 1];
1225  }
1226  else {
1227  rect.boundary = APoints->y[i + 1];
1228  rect.boundary = APoints->y[i];
1229  }
1230  if (APoints->z[i] <= APoints->z[i + 1]) {
1231  rect.boundary = APoints->z[i];
1232  rect.boundary = APoints->z[i + 1];
1233  }
1234  else {
1235  rect.boundary = APoints->z[i + 1];
1236  rect.boundary = APoints->z[i];
1237  }
1238
1239  RTreeSearch(RTree, &rect, (void *)find_cross, &i); /* A segment number from 0 */
1240
1241  if (cross_found) {
1242  break;
1243  }
1244  }
1245
1246  /* Free RTree */
1247  RTreeDestroyNode(RTree);
1248
1249  return cross_found;
1250 }
1251
1265 int
1266 Vect_line_get_intersections(struct line_pnts *APoints,
1267  struct line_pnts *BPoints,
1268  struct line_pnts *IPoints, int with_z)
1269 {
1270  int ret;
1271
1272  IPnts = IPoints;
1273  ret = Vect_line_check_intersection(APoints, BPoints, with_z);
1274
1275  return ret;
1276 }
RectReal boundary[NUMSIDES]
Definition: index.h:42
#define D
int l
double distance
struct line_pnts * Vect_new_line_struct()
Creates and initializes a struct line_pnts.
Definition: line.c:57
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 Node *N, struct Rect *R, SearchHitCallback shcb, void *cbarg)
Definition: index.h:56
void RTreeDestroyNode(struct Node *)
Definition: node.c:217
int Vect_line_check_intersection(struct line_pnts *APoints, struct line_pnts *BPoints, int with_z)
Check if 2 lines intersect.
int RTreeInsertRect(struct Rect *R, int Tid, struct Node **Root, int Level)
int y
Definition: plot.c:34
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:168
#define D2
int dig_line_degenerate(struct line_pnts *points)
Definition: angle.c:180
int Vect_copy_xyz_to_pnts(struct line_pnts *Points, double *x, double *y, double *z, int n)
Copy points from array to line structure.
Definition: line.c:115
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.
#define D1
return NULL
Definition: dbfopen.c:1394
G_warning("category support for [%s] in mapset [%s] %s", name, mapset, type)
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: gis/debug.c:51
int Vect_line_intersection(struct line_pnts *APoints, struct line_pnts *BPoints, struct line_pnts ***ALines, struct line_pnts ***BLines, int *nalines, int *nblines, int with_z)
Intersect 2 lines.
struct Node * RTreeNewIndex(void)
int Vect_line_distance(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 line distance.
Definition: line.c:631
Definition: index.h:40
int Vect_destroy_line_struct(struct line_pnts *p)
Frees all memory associated with a struct line_pnts, including the struct itself. ...
Definition: line.c:90