GRASS 8 Programmer's Manual 8.6.0dev(2026)-ddeab64dbf
Loading...
Searching...
No Matches
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
34mpf_t p11, p12, p21, p22, t1, t2;
35mpf_t dd, dda, ddb, delta;
36mpf_t rra, rrb;
37
38int initialized = 0;
39
40void 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 Calculates:
65 |a11-b11 a12-b12|
66 |a21-b21 a22-b22|
67 */
68void 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
91void 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 */
101int 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 */
361int 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
692int 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 */
928int 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
949int 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
AMI_STREAM
Definition ami_stream.h:153
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:153
MAX
#define MAX(a, b)
Definition gis.h:148
b
double b
Definition r_raster.c:39
t
double t
Definition r_raster.c:39