GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-cf81aadc39
e_intersect.c
Go to the documentation of this file.
1 /*!
2  \file lib/vector/Vlib/e_intersect.c
3 
4  \brief Vector library - intersection (lower level functions)
5 
6  Higher level functions for reading/writing/manipulating vectors.
7 
8  (C) 2008-2009 by the GRASS Development Team
9 
10  This program is free software under the GNU General Public License
11  (>=v2). Read the file COPYING that comes with GRASS for details.
12 
13  \author Rewritten by Rosen Matev (Google Summer of Code 2008)
14  */
15 
16 #include <stdlib.h>
17 #include <math.h>
18 
19 #include <grass/gis.h>
20 
21 #include "e_intersect.h"
22 
23 #define SWAP(a, b) \
24  { \
25  double t = a; \
26  a = b; \
27  b = t; \
28  }
29 #define D ((ax2 - ax1) * (by1 - by2) - (ay2 - ay1) * (bx1 - bx2))
30 #define DA ((bx1 - ax1) * (by1 - by2) - (by1 - ay1) * (bx1 - bx2))
31 #define DB ((ax2 - ax1) * (by1 - ay1) - (ay2 - ay1) * (bx1 - ax1))
32 
33 #ifdef ASDASDASFDSAFFDAS
34 mpf_t p11, p12, p21, p22, t1, t2;
35 mpf_t dd, dda, ddb, delta;
36 mpf_t rra, rrb;
37 
38 int initialized = 0;
39 
40 void initialize_mpf_vars()
41 {
42  mpf_set_default_prec(2048);
43 
44  mpf_init(p11);
45  mpf_init(p12);
46  mpf_init(p21);
47  mpf_init(p22);
48 
49  mpf_init(t1);
50  mpf_init(t2);
51 
52  mpf_init(dd);
53  mpf_init(dda);
54  mpf_init(ddb);
55  mpf_init(delta);
56 
57  mpf_init(rra);
58  mpf_init(rrb);
59 
60  initialized = 1;
61 }
62 
63 /*
64  Caclulates:
65  |a11-b11 a12-b12|
66  |a21-b21 a22-b22|
67  */
68 void det22(mpf_t rop, double a11, double b11, double a12, double b12,
69  double a21, double b21, double a22, double b22)
70 {
71  mpf_set_d(t1, a11);
72  mpf_set_d(t2, b11);
73  mpf_sub(p11, t1, t2);
74  mpf_set_d(t1, a12);
75  mpf_set_d(t2, b12);
76  mpf_sub(p12, t1, t2);
77  mpf_set_d(t1, a21);
78  mpf_set_d(t2, b21);
79  mpf_sub(p21, t1, t2);
80  mpf_set_d(t1, a22);
81  mpf_set_d(t2, b22);
82  mpf_sub(p22, t1, t2);
83 
84  mpf_mul(t1, p11, p22);
85  mpf_mul(t2, p12, p21);
86  mpf_sub(rop, t1, t2);
87 
88  return;
89 }
90 
91 void swap(double *a, double *b)
92 {
93  double t = *a;
94 
95  *a = *b;
96  *b = t;
97  return;
98 }
99 
100 /* multi-precision version */
101 int segment_intersection_2d_e(double ax1, double ay1, double ax2, double ay2,
102  double bx1, double by1, double bx2, double by2,
103  double *x1, double *y1, double *x2, double *y2)
104 {
105  double t;
106 
107  double max_ax, min_ax, max_ay, min_ay;
108 
109  double max_bx, min_bx, max_by, min_by;
110 
111  int sgn_d, sgn_da, sgn_db;
112 
113  int vertical;
114 
115  int f11, f12, f21, f22;
116 
117  mp_exp_t exp;
118 
119  char *s;
120 
121  if (!initialized)
122  initialize_mpf_vars();
123 
124  /* TODO: Works for points ? */
125  G_debug(3, "segment_intersection_2d_e()");
126  G_debug(4, " ax1 = %.18f, ay1 = %.18f", ax1, ay1);
127  G_debug(4, " ax2 = %.18f, ay2 = %.18f", ax2, ay2);
128  G_debug(4, " bx1 = %.18f, by1 = %.18f", bx1, by1);
129  G_debug(4, " bx2 = %.18f, by2 = %.18f", bx2, by2);
130 
131  f11 = ((ax1 == bx1) && (ay1 == by1));
132  f12 = ((ax1 == bx2) && (ay1 == by2));
133  f21 = ((ax2 == bx1) && (ay2 == by1));
134  f22 = ((ax2 == bx2) && (ay2 == by2));
135 
136  /* Check for identical segments */
137  if ((f11 && f22) || (f12 && f21)) {
138  G_debug(3, " identical segments");
139  *x1 = ax1;
140  *y1 = ay1;
141  *x2 = ax2;
142  *y2 = ay2;
143  return 5;
144  }
145  /* Check for identical endpoints */
146  if (f11 || f12) {
147  G_debug(3, " connected by endpoints");
148  *x1 = ax1;
149  *y1 = ay1;
150  return 1;
151  }
152  if (f21 || f22) {
153  G_debug(3, " connected by endpoints");
154  *x1 = ax2;
155  *y1 = ay2;
156  return 1;
157  }
158 
159  if ((MAX(ax1, ax2) < MIN(bx1, bx2)) || (MAX(bx1, bx2) < MIN(ax1, ax2))) {
160  G_debug(3, " no intersection (disjoint bounding boxes)");
161  return 0;
162  }
163  if ((MAX(ay1, ay2) < MIN(by1, by2)) || (MAX(by1, by2) < MIN(ay1, ay2))) {
164  G_debug(3, " no intersection (disjoint bounding boxes)");
165  return 0;
166  }
167 
168  det22(dd, ax2, ax1, bx1, bx2, ay2, ay1, by1, by2);
169  sgn_d = mpf_sgn(dd);
170  if (sgn_d != 0) {
171  G_debug(3, " general position");
172 
173  det22(dda, bx1, ax1, bx1, bx2, by1, ay1, by1, by2);
174  sgn_da = mpf_sgn(dda);
175 
176  /*mpf_div(rra, dda, dd);
177  mpf_div(rrb, ddb, dd);
178  s = mpf_get_str(NULL, &exp, 10, 40, rra);
179  G_debug(4, " ra = %sE%d", (s[0]==0)?"0":s, exp);
180  s = mpf_get_str(NULL, &exp, 10, 24, rrb);
181  G_debug(4, " rb = %sE%d", (s[0]==0)?"0":s, exp);
182  */
183 
184  if (sgn_d > 0) {
185  if ((sgn_da < 0) || (mpf_cmp(dda, dd) > 0)) {
186  G_debug(3, " no intersection");
187  return 0;
188  }
189 
190  det22(ddb, ax2, ax1, bx1, ax1, ay2, ay1, by1, ay1);
191  sgn_db = mpf_sgn(ddb);
192  if ((sgn_db < 0) || (mpf_cmp(ddb, dd) > 0)) {
193  G_debug(3, " no intersection");
194  return 0;
195  }
196  }
197  else { /* if sgn_d < 0 */
198  if ((sgn_da > 0) || (mpf_cmp(dda, dd) < 0)) {
199  G_debug(3, " no intersection");
200  return 0;
201  }
202 
203  det22(ddb, ax2, ax1, bx1, ax1, ay2, ay1, by1, ay1);
204  sgn_db = mpf_sgn(ddb);
205  if ((sgn_db > 0) || (mpf_cmp(ddb, dd) < 0)) {
206  G_debug(3, " no intersection");
207  return 0;
208  }
209  }
210 
211  /*G_debug(3, " ra=%.17g rb=%.17g", mpf_get_d(dda)/mpf_get_d(dd),
212  * mpf_get_d(ddb)/mpf_get_d(dd)); */
213  /*G_debug(3, " sgn_d=%d sgn_da=%d sgn_db=%d cmp(dda,dd)=%d
214  * cmp(ddb,dd)=%d", sgn_d, sgn_da, sgn_db, mpf_cmp(dda, dd),
215  * mpf_cmp(ddb, dd)); */
216 
217  mpf_set_d(delta, ax2 - ax1);
218  mpf_mul(t1, dda, delta);
219  mpf_div(t2, t1, dd);
220  *x1 = ax1 + mpf_get_d(t2);
221 
222  mpf_set_d(delta, ay2 - ay1);
223  mpf_mul(t1, dda, delta);
224  mpf_div(t2, t1, dd);
225  *y1 = ay1 + mpf_get_d(t2);
226 
227  G_debug(3, " intersection %.16g, %.16g", *x1, *y1);
228  return 1;
229  }
230 
231  /* segments are parallel or collinear */
232  det22(dda, bx1, ax1, bx1, bx2, by1, ay1, by1, by2);
233  sgn_da = mpf_sgn(dda);
234  if (sgn_da != 0) {
235  /* segments are parallel */
236  G_debug(3, " parallel segments");
237  return 0;
238  }
239 
240  /* segments are colinear. check for overlap */
241 
242  /* swap endpoints if needed */
243  /* if segments are vertical, we swap x-coords with y-coords */
244  vertical = 0;
245  if (ax1 > ax2) {
246  SWAP(ax1, ax2);
247  SWAP(ay1, ay2);
248  }
249  else if (ax1 == ax2) {
250  vertical = 1;
251  if (ay1 > ay2)
252  SWAP(ay1, ay2);
253  SWAP(ax1, ay1);
254  SWAP(ax2, ay2);
255  }
256  if (bx1 > bx2) {
257  SWAP(bx1, bx2);
258  SWAP(by1, by2);
259  }
260  else if (bx1 == bx2) {
261  if (by1 > by2)
262  SWAP(by1, by2);
263  SWAP(bx1, by1);
264  SWAP(bx2, by2);
265  }
266 
267  G_debug(3, " collinear segments");
268 
269  if ((bx2 < ax1) || (bx1 > ax2)) {
270  G_debug(3, " no intersection");
271  return 0;
272  }
273 
274  /* there is overlap or connected end points */
275  G_debug(3, " overlap");
276 
277  /* a contains b */
278  if ((ax1 < bx1) && (ax2 > bx2)) {
279  G_debug(3, " a contains b");
280  if (!vertical) {
281  *x1 = bx1;
282  *y1 = by1;
283  *x2 = bx2;
284  *y2 = by2;
285  }
286  else {
287  *x1 = by1;
288  *y1 = bx1;
289  *x2 = by2;
290  *y2 = bx2;
291  }
292  return 3;
293  }
294 
295  /* b contains a */
296  if ((ax1 > bx1) && (ax2 < bx2)) {
297  G_debug(3, " b contains a");
298  if (!vertical) {
299  *x1 = bx1;
300  *y1 = by1;
301  *x2 = bx2;
302  *y2 = by2;
303  }
304  else {
305  *x1 = by1;
306  *y1 = bx1;
307  *x2 = by2;
308  *y2 = bx2;
309  }
310  return 4;
311  }
312 
313  /* general overlap, 2 intersection points */
314  G_debug(3, " partial overlap");
315  if ((bx1 > ax1) && (bx1 < ax2)) { /* b1 is in a */
316  if (!vertical) {
317  *x1 = bx1;
318  *y1 = by1;
319  *x2 = ax2;
320  *y2 = ay2;
321  }
322  else {
323  *x1 = by1;
324  *y1 = bx1;
325  *x2 = ay2;
326  *y2 = ax2;
327  }
328  return 2;
329  }
330  if ((bx2 > ax1) && (bx2 < ax2)) { /* b2 is in a */
331  if (!vertical) {
332  *x1 = bx2;
333  *y1 = by2;
334  *x2 = ax1;
335  *y2 = ay1;
336  }
337  else {
338  *x1 = by2;
339  *y1 = bx2;
340  *x2 = ay1;
341  *y2 = ax1;
342  }
343  return 2;
344  }
345 
346  /* should not be reached */
347  G_warning(("segment_intersection_2d() ERROR (should not be reached)"));
348  G_warning("%.16g %.16g", ax1, ay1);
349  G_warning("%.16g %.16g", ax2, ay2);
350  G_warning("x");
351  G_warning("%.16g %.16g", bx1, by1);
352  G_warning("%.16g %.16g", bx2, by2);
353 
354  return 0;
355 }
356 #endif
357 
358 /* OLD */
359 /* tolerance aware version */
360 /* TODO: fix all ==s left */
361 int segment_intersection_2d_tol(double ax1, double ay1, double ax2, double ay2,
362  double bx1, double by1, double bx2, double by2,
363  double *x1, double *y1, double *x2, double *y2,
364  double tol)
365 {
366  double tola, tolb;
367 
368  double d, d1, d2, ra, rb, t;
369 
370  int switched = 0;
371 
372  /* TODO: Works for points ? */
373  G_debug(4, "segment_intersection_2d()");
374  G_debug(4, " ax1 = %.18f, ay1 = %.18f", ax1, ay1);
375  G_debug(4, " ax2 = %.18f, ay2 = %.18f", ax2, ay2);
376  G_debug(4, " bx1 = %.18f, by1 = %.18f", bx1, by1);
377  G_debug(4, " bx2 = %.18f, by2 = %.18f", bx2, by2);
378 
379  /* Check identical lines */
380  if ((FEQUAL(ax1, bx1, tol) && FEQUAL(ay1, by1, tol) &&
381  FEQUAL(ax2, bx2, tol) && FEQUAL(ay2, by2, tol)) ||
382  (FEQUAL(ax1, bx2, tol) && FEQUAL(ay1, by2, tol) &&
383  FEQUAL(ax2, bx1, tol) && FEQUAL(ay2, by1, tol))) {
384  G_debug(2, " -> identical segments");
385  *x1 = ax1;
386  *y1 = ay1;
387  *x2 = ax2;
388  *y2 = ay2;
389  return 5;
390  }
391 
392  /* 'Sort' lines by x1, x2, y1, y2 */
393  if (bx1 < ax1)
394  switched = 1;
395  else if (bx1 == ax1) {
396  if (bx2 < ax2)
397  switched = 1;
398  else if (bx2 == ax2) {
399  if (by1 < ay1)
400  switched = 1;
401  else if (by1 == ay1) {
402  if (by2 < ay2)
403  switched = 1; /* by2 != ay2 (would be identical */
404  }
405  }
406  }
407 
408  if (switched) {
409  t = ax1;
410  ax1 = bx1;
411  bx1 = t;
412  t = ay1;
413  ay1 = by1;
414  by1 = t;
415  t = ax2;
416  ax2 = bx2;
417  bx2 = t;
418  t = ay2;
419  ay2 = by2;
420  by2 = t;
421  }
422 
423  d = (ax2 - ax1) * (by1 - by2) - (ay2 - ay1) * (bx1 - bx2);
424  d1 = (bx1 - ax1) * (by1 - by2) - (by1 - ay1) * (bx1 - bx2);
425  d2 = (ax2 - ax1) * (by1 - ay1) - (ay2 - ay1) * (bx1 - ax1);
426 
427  G_debug(2, " d = %.18g", d);
428  G_debug(2, " d1 = %.18g", d1);
429  G_debug(2, " d2 = %.18g", d2);
430 
431  tola = tol / MAX(fabs(ax2 - ax1), fabs(ay2 - ay1));
432  tolb = tol / MAX(fabs(bx2 - bx1), fabs(by2 - by1));
433  G_debug(2, " tol = %.18g", tol);
434  G_debug(2, " tola = %.18g", tola);
435  G_debug(2, " tolb = %.18g", tolb);
436  if (!FZERO(d, tol)) {
437  ra = d1 / d;
438  rb = d2 / d;
439 
440  G_debug(2, " not parallel/collinear: ra = %.18g", ra);
441  G_debug(2, " rb = %.18g", rb);
442 
443  if ((ra <= -tola) || (ra >= 1 + tola) || (rb <= -tolb) ||
444  (rb >= 1 + tolb)) {
445  G_debug(2, " no intersection");
446  return 0;
447  }
448 
449  ra = MIN(MAX(ra, 0), 1);
450  *x1 = ax1 + ra * (ax2 - ax1);
451  *y1 = ay1 + ra * (ay2 - ay1);
452 
453  G_debug(2, " intersection %.18f, %.18f", *x1, *y1);
454  return 1;
455  }
456 
457  /* segments are parallel or collinear */
458  G_debug(3, " -> parallel/collinear");
459 
460  if ((!FZERO(d1, tol)) || (!FZERO(d2, tol))) { /* lines are parallel */
461  G_debug(2, " -> parallel");
462  return 0;
463  }
464 
465  /* segments are colinear. check for overlap */
466 
467  /* aa = adx*adx + ady*ady;
468  bb = bdx*bdx + bdy*bdy;
469 
470  t = (ax1-bx1)*bdx + (ay1-by1)*bdy; */
471 
472  /* Collinear vertical */
473  /* original code assumed lines were not both vertical
474  * so there is a special case if they are */
475  if (FEQUAL(ax1, ax2, tol) && FEQUAL(bx1, bx2, tol) &&
476  FEQUAL(ax1, bx1, tol)) {
477  G_debug(2, " -> collinear vertical");
478  if (ay1 > ay2) {
479  t = ay1;
480  ay1 = ay2;
481  ay2 = t;
482  } /* to be sure that ay1 < ay2 */
483  if (by1 > by2) {
484  t = by1;
485  by1 = by2;
486  by2 = t;
487  } /* to be sure that by1 < by2 */
488  if (ay1 > by2 || ay2 < by1) {
489  G_debug(2, " -> no intersection");
490  return 0;
491  }
492 
493  /* end points */
494  if (FEQUAL(ay1, by2, tol)) {
495  *x1 = ax1;
496  *y1 = ay1;
497  G_debug(2, " -> connected by end points");
498  return 1; /* endpoints only */
499  }
500  if (FEQUAL(ay2, by1, tol)) {
501  *x1 = ax2;
502  *y1 = ay2;
503  G_debug(2, " -> connected by end points");
504  return 1; /* endpoints only */
505  }
506 
507  /* general overlap */
508  G_debug(3, " -> vertical overlap");
509  /* a contains b */
510  if (ay1 <= by1 && ay2 >= by2) {
511  G_debug(2, " -> a contains b");
512  *x1 = bx1;
513  *y1 = by1;
514  *x2 = bx2;
515  *y2 = by2;
516  if (!switched)
517  return 3;
518  else
519  return 4;
520  }
521  /* b contains a */
522  if (ay1 >= by1 && ay2 <= by2) {
523  G_debug(2, " -> b contains a");
524  *x1 = ax1;
525  *y1 = ay1;
526  *x2 = ax2;
527  *y2 = ay2;
528  if (!switched)
529  return 4;
530  else
531  return 3;
532  }
533 
534  /* general overlap, 2 intersection points */
535  G_debug(2, " -> partial overlap");
536  if (by1 > ay1 && by1 < ay2) { /* b1 in a */
537  if (!switched) {
538  *x1 = bx1;
539  *y1 = by1;
540  *x2 = ax2;
541  *y2 = ay2;
542  }
543  else {
544  *x1 = ax2;
545  *y1 = ay2;
546  *x2 = bx1;
547  *y2 = by1;
548  }
549  return 2;
550  }
551 
552  if (by2 > ay1 && by2 < ay2) { /* b2 in a */
553  if (!switched) {
554  *x1 = bx2;
555  *y1 = by2;
556  *x2 = ax1;
557  *y2 = ay1;
558  }
559  else {
560  *x1 = ax1;
561  *y1 = ay1;
562  *x2 = bx2;
563  *y2 = by2;
564  }
565  return 2;
566  }
567 
568  /* should not be reached */
569  G_warning((
570  "Vect_segment_intersection() ERROR (collinear vertical segments)"));
571  G_warning("%.15g %.15g", ax1, ay1);
572  G_warning("%.15g %.15g", ax2, ay2);
573  G_warning("x");
574  G_warning("%.15g %.15g", bx1, by1);
575  G_warning("%.15g %.15g", bx2, by2);
576  return 0;
577  }
578 
579  G_debug(2, " -> collinear non vertical");
580 
581  /* Collinear non vertical */
582  if ((bx1 > ax1 && bx2 > ax1 && bx1 > ax2 && bx2 > ax2) ||
583  (bx1 < ax1 && bx2 < ax1 && bx1 < ax2 && bx2 < ax2)) {
584  G_debug(2, " -> no intersection");
585  return 0;
586  }
587 
588  /* there is overlap or connected end points */
589  G_debug(2, " -> overlap/connected end points");
590 
591  /* end points */
592  if ((ax1 == bx1 && ay1 == by1) || (ax1 == bx2 && ay1 == by2)) {
593  *x1 = ax1;
594  *y1 = ay1;
595  G_debug(2, " -> connected by end points");
596  return 1;
597  }
598  if ((ax2 == bx1 && ay2 == by1) || (ax2 == bx2 && ay2 == by2)) {
599  *x1 = ax2;
600  *y1 = ay2;
601  G_debug(2, " -> connected by end points");
602  return 1;
603  }
604 
605  if (ax1 > ax2) {
606  t = ax1;
607  ax1 = ax2;
608  ax2 = t;
609  t = ay1;
610  ay1 = ay2;
611  ay2 = t;
612  } /* to be sure that ax1 < ax2 */
613  if (bx1 > bx2) {
614  t = bx1;
615  bx1 = bx2;
616  bx2 = t;
617  t = by1;
618  by1 = by2;
619  by2 = t;
620  } /* to be sure that bx1 < bx2 */
621 
622  /* a contains b */
623  if (ax1 <= bx1 && ax2 >= bx2) {
624  G_debug(2, " -> a contains b");
625  *x1 = bx1;
626  *y1 = by1;
627  *x2 = bx2;
628  *y2 = by2;
629  if (!switched)
630  return 3;
631  else
632  return 4;
633  }
634  /* b contains a */
635  if (ax1 >= bx1 && ax2 <= bx2) {
636  G_debug(2, " -> b contains a");
637  *x1 = ax1;
638  *y1 = ay1;
639  *x2 = ax2;
640  *y2 = ay2;
641  if (!switched)
642  return 4;
643  else
644  return 3;
645  }
646 
647  /* general overlap, 2 intersection points (lines are not vertical) */
648  G_debug(2, " -> partial overlap");
649  if (bx1 > ax1 && bx1 < ax2) { /* b1 is in a */
650  if (!switched) {
651  *x1 = bx1;
652  *y1 = by1;
653  *x2 = ax2;
654  *y2 = ay2;
655  }
656  else {
657  *x1 = ax2;
658  *y1 = ay2;
659  *x2 = bx1;
660  *y2 = by1;
661  }
662  return 2;
663  }
664  if (bx2 > ax1 && bx2 < ax2) { /* b2 is in a */
665  if (!switched) {
666  *x1 = bx2;
667  *y1 = by2;
668  *x2 = ax1;
669  *y2 = ay1;
670  }
671  else {
672  *x1 = ax1;
673  *y1 = ay1;
674  *x2 = bx2;
675  *y2 = by2;
676  }
677  return 2;
678  }
679 
680  /* should not be reached */
681  G_warning(
682  ("segment_intersection_2d() ERROR (collinear non vertical segments)"));
683  G_warning("%.15g %.15g", ax1, ay1);
684  G_warning("%.15g %.15g", ax2, ay2);
685  G_warning("x");
686  G_warning("%.15g %.15g", bx1, by1);
687  G_warning("%.15g %.15g", bx2, by2);
688 
689  return 0;
690 }
691 
692 int segment_intersection_2d(double ax1, double ay1, double ax2, double ay2,
693  double bx1, double by1, double bx2, double by2,
694  double *x1, double *y1, double *x2, double *y2)
695 {
696  const int DLEVEL = 4;
697 
698  int vertical;
699 
700  int f11, f12, f21, f22;
701 
702  double d, da, db;
703 
704  /* TODO: Works for points ? */
705  G_debug(DLEVEL, "segment_intersection_2d()");
706  G_debug(4, " ax1 = %.18f, ay1 = %.18f", ax1, ay1);
707  G_debug(4, " ax2 = %.18f, ay2 = %.18f", ax2, ay2);
708  G_debug(4, " bx1 = %.18f, by1 = %.18f", bx1, by1);
709  G_debug(4, " bx2 = %.18f, by2 = %.18f", bx2, by2);
710 
711  f11 = ((ax1 == bx1) && (ay1 == by1));
712  f12 = ((ax1 == bx2) && (ay1 == by2));
713  f21 = ((ax2 == bx1) && (ay2 == by1));
714  f22 = ((ax2 == bx2) && (ay2 == by2));
715 
716  /* Check for identical segments */
717  if ((f11 && f22) || (f12 && f21)) {
718  G_debug(DLEVEL, " identical segments");
719  *x1 = ax1;
720  *y1 = ay1;
721  *x2 = ax2;
722  *y2 = ay2;
723  return 5;
724  }
725  /* Check for identical endpoints */
726  if (f11 || f12) {
727  G_debug(DLEVEL, " connected by endpoints");
728  *x1 = ax1;
729  *y1 = ay1;
730  return 1;
731  }
732  if (f21 || f22) {
733  G_debug(DLEVEL, " connected by endpoints");
734  *x1 = ax2;
735  *y1 = ay2;
736  return 1;
737  }
738 
739  if ((MAX(ax1, ax2) < MIN(bx1, bx2)) || (MAX(bx1, bx2) < MIN(ax1, ax2))) {
740  G_debug(DLEVEL, " no intersection (disjoint bounding boxes)");
741  return 0;
742  }
743  if ((MAX(ay1, ay2) < MIN(by1, by2)) || (MAX(by1, by2) < MIN(ay1, ay2))) {
744  G_debug(DLEVEL, " no intersection (disjoint bounding boxes)");
745  return 0;
746  }
747 
748  /* swap endpoints if needed */
749  /* if segments are vertical, we swap x-coords with y-coords */
750  vertical = 0;
751  if (ax1 > ax2) {
752  SWAP(ax1, ax2);
753  SWAP(ay1, ay2);
754  }
755  else if (ax1 == ax2) {
756  vertical = 1;
757  if (ay1 > ay2)
758  SWAP(ay1, ay2);
759  SWAP(ax1, ay1);
760  SWAP(ax2, ay2);
761  }
762  if (bx1 > bx2) {
763  SWAP(bx1, bx2);
764  SWAP(by1, by2);
765  }
766  else if (bx1 == bx2) {
767  if (by1 > by2)
768  SWAP(by1, by2);
769  SWAP(bx1, by1);
770  SWAP(bx2, by2);
771  }
772 
773  d = D;
774  if (d != 0) {
775  G_debug(DLEVEL, " general position");
776 
777  da = DA;
778 
779  /*mpf_div(rra, dda, dd);
780  mpf_div(rrb, ddb, dd);
781  s = mpf_get_str(NULL, &exp, 10, 40, rra);
782  G_debug(4, " ra = %sE%d", (s[0]==0)?"0":s, exp);
783  s = mpf_get_str(NULL, &exp, 10, 24, rrb);
784  G_debug(4, " rb = %sE%d", (s[0]==0)?"0":s, exp);
785  */
786 
787  if (d > 0) {
788  if ((da < 0) || (da > d)) {
789  G_debug(DLEVEL, " no intersection");
790  return 0;
791  }
792 
793  db = DB;
794  if ((db < 0) || (db > d)) {
795  G_debug(DLEVEL, " no intersection");
796  return 0;
797  }
798  }
799  else { /* if d < 0 */
800  if ((da > 0) || (da < d)) {
801  G_debug(DLEVEL, " no intersection");
802  return 0;
803  }
804 
805  db = DB;
806  if ((db > 0) || (db < d)) {
807  G_debug(DLEVEL, " no intersection");
808  return 0;
809  }
810  }
811 
812  /*G_debug(DLEVEL, " ra=%.17g rb=%.17g",
813  * mpf_get_d(dda)/mpf_get_d(dd), mpf_get_d(ddb)/mpf_get_d(dd)); */
814  /*G_debug(DLEVEL, " sgn_d=%d sgn_da=%d sgn_db=%d cmp(dda,dd)=%d
815  * cmp(ddb,dd)=%d", sgn_d, sgn_da, sgn_db, mpf_cmp(dda, dd),
816  * mpf_cmp(ddb, dd)); */
817 
818  *x1 = ax1 + (ax2 - ax1) * da / d;
819  *y1 = ay1 + (ay2 - ay1) * da / d;
820 
821  G_debug(DLEVEL, " intersection %.16g, %.16g", *x1, *y1);
822  return 1;
823  }
824 
825  /* segments are parallel or collinear */
826  da = DA;
827  db = DB;
828  if ((da != 0) || (db != 0)) {
829  /* segments are parallel */
830  G_debug(DLEVEL, " parallel segments");
831  return 0;
832  }
833 
834  /* segments are colinear. check for overlap */
835 
836  G_debug(DLEVEL, " collinear segments");
837 
838  if ((bx2 < ax1) || (bx1 > ax2)) {
839  G_debug(DLEVEL, " no intersection");
840  return 0;
841  }
842 
843  /* there is overlap or connected end points */
844  G_debug(DLEVEL, " overlap");
845 
846  /* a contains b */
847  if ((ax1 < bx1) && (ax2 > bx2)) {
848  G_debug(DLEVEL, " a contains b");
849  if (!vertical) {
850  *x1 = bx1;
851  *y1 = by1;
852  *x2 = bx2;
853  *y2 = by2;
854  }
855  else {
856  *x1 = by1;
857  *y1 = bx1;
858  *x2 = by2;
859  *y2 = bx2;
860  }
861  return 3;
862  }
863 
864  /* b contains a */
865  if ((ax1 > bx1) && (ax2 < bx2)) {
866  G_debug(DLEVEL, " b contains a");
867  if (!vertical) {
868  *x1 = bx1;
869  *y1 = by1;
870  *x2 = bx2;
871  *y2 = by2;
872  }
873  else {
874  *x1 = by1;
875  *y1 = bx1;
876  *x2 = by2;
877  *y2 = bx2;
878  }
879  return 4;
880  }
881 
882  /* general overlap, 2 intersection points */
883  G_debug(DLEVEL, " partial overlap");
884  if ((bx1 > ax1) && (bx1 < ax2)) { /* b1 is in a */
885  if (!vertical) {
886  *x1 = bx1;
887  *y1 = by1;
888  *x2 = ax2;
889  *y2 = ay2;
890  }
891  else {
892  *x1 = by1;
893  *y1 = bx1;
894  *x2 = ay2;
895  *y2 = ax2;
896  }
897  return 2;
898  }
899  if ((bx2 > ax1) && (bx2 < ax2)) { /* b2 is in a */
900  if (!vertical) {
901  *x1 = bx2;
902  *y1 = by2;
903  *x2 = ax1;
904  *y2 = ay1;
905  }
906  else {
907  *x1 = by2;
908  *y1 = bx2;
909  *x2 = ay1;
910  *y2 = ax1;
911  }
912  return 2;
913  }
914 
915  /* should not be reached */
916  G_warning(("segment_intersection_2d() ERROR (should not be reached)"));
917  G_warning("%.16g %.16g", ax1, ay1);
918  G_warning("%.16g %.16g", ax2, ay2);
919  G_warning("x");
920  G_warning("%.16g %.16g", bx1, by1);
921  G_warning("%.16g %.16g", bx2, by2);
922 
923  return 0;
924 }
925 
926 #define N 52 /* double's mantisa size in bits */
927 /* a and b are different in at most <bits> significant digits */
928 int almost_equal(double a, double b, int bits)
929 {
930  int ea, eb, e;
931 
932  if (a == b)
933  return 1;
934 
935  if (a == 0 || b == 0) {
936  /* return (0 < -N+bits); */
937  return (bits > N);
938  }
939 
940  frexp(a, &ea);
941  frexp(b, &eb);
942  if (ea != eb)
943  return (bits > N + abs(ea - eb));
944  frexp(a - b, &e);
945  return (e < ea - N + bits);
946 }
947 
948 #ifdef ASDASDFASFEAS
949 int segment_intersection_2d_test(double ax1, double ay1, double ax2, double ay2,
950  double bx1, double by1, double bx2, double by2,
951  double *x1, double *y1, double *x2, double *y2)
952 {
953  double t;
954 
955  double max_ax, min_ax, max_ay, min_ay;
956 
957  double max_bx, min_bx, max_by, min_by;
958 
959  int sgn_d, sgn_da, sgn_db;
960 
961  int vertical;
962 
963  int f11, f12, f21, f22;
964 
965  mp_exp_t exp;
966 
967  char *s;
968 
969  double d, da, db, ra, rb;
970 
971  if (!initialized)
972  initialize_mpf_vars();
973 
974  /* TODO: Works for points ? */
975  G_debug(3, "segment_intersection_2d_test()");
976  G_debug(3, " ax1 = %.18e, ay1 = %.18e", ax1, ay1);
977  G_debug(3, " ax2 = %.18e, ay2 = %.18e", ax2, ay2);
978  G_debug(3, " bx1 = %.18e, by1 = %.18e", bx1, by1);
979  G_debug(3, " bx2 = %.18e, by2 = %.18e", bx2, by2);
980 
981  f11 = ((ax1 == bx1) && (ay1 == by1));
982  f12 = ((ax1 == bx2) && (ay1 == by2));
983  f21 = ((ax2 == bx1) && (ay2 == by1));
984  f22 = ((ax2 == bx2) && (ay2 == by2));
985 
986  /* Check for identical segments */
987  if ((f11 && f22) || (f12 && f21)) {
988  G_debug(4, " identical segments");
989  *x1 = ax1;
990  *y1 = ay1;
991  *x2 = ax2;
992  *y2 = ay2;
993  return 5;
994  }
995  /* Check for identical endpoints */
996  if (f11 || f12) {
997  G_debug(4, " connected by endpoints");
998  *x1 = ax1;
999  *y1 = ay1;
1000  return 1;
1001  }
1002  if (f21 || f22) {
1003  G_debug(4, " connected by endpoints");
1004  *x1 = ax2;
1005  *y1 = ay2;
1006  return 1;
1007  }
1008 
1009  if ((MAX(ax1, ax2) < MIN(bx1, bx2)) || (MAX(bx1, bx2) < MIN(ax1, ax2))) {
1010  G_debug(4, " no intersection (disjoint bounding boxes)");
1011  return 0;
1012  }
1013  if ((MAX(ay1, ay2) < MIN(by1, by2)) || (MAX(by1, by2) < MIN(ay1, ay2))) {
1014  G_debug(4, " no intersection (disjoint bounding boxes)");
1015  return 0;
1016  }
1017 
1018  d = (ax2 - ax1) * (by1 - by2) - (ay2 - ay1) * (bx1 - bx2);
1019  da = (bx1 - ax1) * (by1 - by2) - (by1 - ay1) * (bx1 - bx2);
1020  db = (ax2 - ax1) * (by1 - ay1) - (ay2 - ay1) * (bx1 - ax1);
1021 
1022  det22(dd, ax2, ax1, bx1, bx2, ay2, ay1, by1, by2);
1023  sgn_d = mpf_sgn(dd);
1024  s = mpf_get_str(NULL, &exp, 10, 40, dd);
1025  G_debug(3, " dd = %sE%d", (s[0] == 0) ? "0" : s, exp);
1026  G_debug(3, " d = %.18E", d);
1027 
1028  if (sgn_d != 0) {
1029  G_debug(3, " general position");
1030 
1031  det22(dda, bx1, ax1, bx1, bx2, by1, ay1, by1, by2);
1032  det22(ddb, ax2, ax1, bx1, ax1, ay2, ay1, by1, ay1);
1033  sgn_da = mpf_sgn(dda);
1034  sgn_db = mpf_sgn(ddb);
1035 
1036  ra = da / d;
1037  rb = db / d;
1038  mpf_div(rra, dda, dd);
1039  mpf_div(rrb, ddb, dd);
1040 
1041  s = mpf_get_str(NULL, &exp, 10, 40, rra);
1042  G_debug(4, " rra = %sE%d", (s[0] == 0) ? "0" : s, exp);
1043  G_debug(4, " ra = %.18E", ra);
1044  s = mpf_get_str(NULL, &exp, 10, 40, rrb);
1045  G_debug(4, " rrb = %sE%d", (s[0] == 0) ? "0" : s, exp);
1046  G_debug(4, " rb = %.18E", rb);
1047 
1048  if (sgn_d > 0) {
1049  if ((sgn_da < 0) || (mpf_cmp(dda, dd) > 0)) {
1050  G_debug(DLEVEL, " no intersection");
1051  return 0;
1052  }
1053 
1054  if ((sgn_db < 0) || (mpf_cmp(ddb, dd) > 0)) {
1055  G_debug(DLEVEL, " no intersection");
1056  return 0;
1057  }
1058  }
1059  else { /* if sgn_d < 0 */
1060  if ((sgn_da > 0) || (mpf_cmp(dda, dd) < 0)) {
1061  G_debug(DLEVEL, " no intersection");
1062  return 0;
1063  }
1064 
1065  if ((sgn_db > 0) || (mpf_cmp(ddb, dd) < 0)) {
1066  G_debug(DLEVEL, " no intersection");
1067  return 0;
1068  }
1069  }
1070 
1071  mpf_set_d(delta, ax2 - ax1);
1072  mpf_mul(t1, dda, delta);
1073  mpf_div(t2, t1, dd);
1074  *x1 = ax1 + mpf_get_d(t2);
1075 
1076  mpf_set_d(delta, ay2 - ay1);
1077  mpf_mul(t1, dda, delta);
1078  mpf_div(t2, t1, dd);
1079  *y1 = ay1 + mpf_get_d(t2);
1080 
1081  G_debug(2, " intersection at:");
1082  G_debug(2, " xx = %.18e", *x1);
1083  G_debug(2, " x = %.18e", ax1 + ra * (ax2 - ax1));
1084  G_debug(2, " yy = %.18e", *y1);
1085  G_debug(2, " y = %.18e", ay1 + ra * (ay2 - ay1));
1086  return 1;
1087  }
1088 
1089  G_debug(3, " parallel/collinear...");
1090  return -1;
1091 }
1092 #endif
NULL
#define NULL
Definition: ccmath.h:32
G_warning
void G_warning(const char *,...) __attribute__((format(printf
G_debug
int G_debug(int, const char *,...) __attribute__((format(printf
N
#define N
Definition: e_intersect.c:926
segment_intersection_2d
int segment_intersection_2d(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2, double *x1, double *y1, double *x2, double *y2)
Definition: e_intersect.c:692
segment_intersection_2d_tol
int segment_intersection_2d_tol(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2, double *x1, double *y1, double *x2, double *y2, double tol)
Definition: e_intersect.c:361
DB
#define DB
Definition: e_intersect.c:31
SWAP
#define SWAP(a, b)
Definition: e_intersect.c:23
DA
#define DA
Definition: e_intersect.c:30
almost_equal
int almost_equal(double a, double b, int bits)
Definition: e_intersect.c:928
D
#define D
Definition: e_intersect.c:29
e_intersect.h
FZERO
#define FZERO(X, TOL)
Definition: e_intersect.h:4
FEQUAL
#define FEQUAL(X, Y, TOL)
Definition: e_intersect.h:5
MIN
#define MIN(a, b)
Definition: gis.h:154
MAX
#define MAX(a, b)
Definition: gis.h:149
b
double b
Definition: r_raster.c:39
t
double t
Definition: r_raster.c:39