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