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