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