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