GRASS GIS 8 Programmer's Manual 8.2.1(2023)-exported
intersect.c
Go to the documentation of this file.
1
2/**************************************************************
3 * find intersection between two lines defined by points on the lines
4 * line segment A is (ax1,ay1) to (ax2,ay2)
5 * line segment B is (bx1,by1) to (bx2,by2)
6 * returns
7 * -1 segment A and B do not intersect (parallel without overlap)
8 * 0 segment A and B do not intersect but extensions do intersect
9 * 1 intersection is a single point
10 * 2 intersection is a line segment (colinear with overlap)
11 * x,y intersection point
12 * ra - ratio that the intersection divides A
13 * rb - ratio that the intersection divides B
14 *
15 * B2
16 * /
17 * /
18 * r=p/(p+q) : A1---p-------*--q------A2
19 * /
20 * /
21 * B1
22 *
23 **************************************************************/
24
25/**************************************************************
26*
27* A point P which lies on line defined by points A1=(x1,y1) and A2=(x2,y2)
28* is given by the equation r * (x2,y2) + (1-r) * (x1,y1).
29* if r is between 0 and 1, p lies between A1 and A2.
30*
31* Suppose points on line (A1, A2) has equation
32* (x,y) = ra * (ax2,ay2) + (1-ra) * (ax1,ay1)
33* or for x and y separately
34* x = ra * ax2 - ra * ax1 + ax1
35* y = ra * ay2 - ra * ay1 + ay1
36* and the points on line (B1, B2) are represented by
37* (x,y) = rb * (bx2,by2) + (1-rb) * (bx1,by1)
38* or for x and y separately
39* x = rb * bx2 - rb * bx1 + bx1
40* y = rb * by2 - rb * by1 + by1
41*
42* when the lines intersect, the point (x,y) has to
43* satisfy a system of 2 equations:
44* ra * ax2 - ra * ax1 + ax1 = rb * bx2 - rb * bx1 + bx1
45* ra * ay2 - ra * ay1 + ay1 = rb * by2 - rb * by1 + by1
46*
47* or
48*
49* (ax2 - ax1) * ra - (bx2 - bx1) * rb = bx1 - ax1
50* (ay2 - ay1) * ra - (by2 - by1) * rb = by1 - ay1
51*
52* by Cramer's method, one can solve this by computing 3
53* determinants of matrices:
54*
55* M = (ax2-ax1) (bx1-bx2)
56* (ay2-ay1) (by1-by2)
57*
58* M1 = (bx1-ax1) (bx1-bx2)
59* (by1-ay1) (by1-by2)
60*
61* M2 = (ax2-ax1) (bx1-ax1)
62* (ay2-ay1) (by1-ay1)
63*
64* Which are exactly the determinants D, D2, D1 below:
65*
66* D ((ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2))
67*
68* D1 ((bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2))
69*
70* D2 ((ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1))
71***********************************************************************/
72
73
74#define D ((ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2))
75#define D1 ((bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2))
76#define D2 ((ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1))
77
78#define SWAP(x,y) {double t; t=x; x=y; y=t;}
79
80int G_intersect_line_segments(double ax1, double ay1, double ax2, double ay2,
81 double bx1, double by1, double bx2, double by2,
82 double *ra, double *rb, double *x, double *y)
83{
84 double d;
85
86 if (ax1 > ax2 || (ax1 == ax2 && ay1 > ay2)) {
87 SWAP(ax1, ax2)
88 SWAP(ay1, ay2)
89 }
90
91 if (bx1 > bx2 || (bx1 == bx2 && by1 > by2)) {
92 SWAP(bx1, bx2)
93 SWAP(by1, by2)
94 }
95
96 d = D;
97
98 if (d) { /* lines are not parallel */
99 *ra = D1 / d;
100 *rb = D2 / d;
101
102 *x = ax1 + (*ra) * (ax2 - ax1);
103 *y = ay1 + (*ra) * (ay2 - ay1);
104 return (*ra >= 0.0 && *ra <= 1.0 && *rb >= 0.0 && *rb <= 1.0);
105 }
106
107 /* lines are parallel */
108 if (D1 || D2) {
109 return -1;
110 }
111
112 /* segments are colinear. check for overlap */
113
114 /* Collinear vertical */
115 if (ax1 == ax2) {
116 if (ay1 > by2) {
117 *x = ax1;
118 *y = ay1;
119 return 0; /* extensions overlap */
120 }
121 if (ay2 < by1) {
122 *x = ax2;
123 *y = ay2;
124 return 0; /* extensions overlap */
125 }
126
127 /* there is overlap */
128
129 if (ay1 == by2) {
130 *x = ax1;
131 *y = ay1;
132 return 1; /* endpoints only */
133 }
134 if (ay2 == by1) {
135 *x = ax2;
136 *y = ay2;
137 return 1; /* endpoints only */
138 }
139
140 /* overlap, no single intersection point */
141 if (ay1 > by1 && ay1 < by2) {
142 *x = ax1;
143 *y = ay1;
144 }
145 else {
146 *x = ax2;
147 *y = ay2;
148 }
149 return 2;
150 }
151 else {
152 if (ax1 > bx2) {
153 *x = ax1;
154 *y = ay1;
155 return 0; /* extensions overlap */
156 }
157 if (ax2 < bx1) {
158 *x = ax2;
159 *y = ay2;
160 return 0; /* extensions overlap */
161 }
162
163 /* there is overlap */
164
165 if (ax1 == bx2) {
166 *x = ax1;
167 *y = ay1;
168 return 1; /* endpoints only */
169 }
170 if (ax2 == bx1) {
171 *x = ax2;
172 *y = ay2;
173 return 1; /* endpoints only */
174 }
175
176 /* overlap, no single intersection point */
177 if (ax1 > bx1 && ax1 < bx2) {
178 *x = ax1;
179 *y = ay1;
180 }
181 else {
182 *x = ax2;
183 *y = ay2;
184 }
185 return 2;
186 }
187
188 return 0; /* should not be reached */
189}
int G_intersect_line_segments(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2, double *ra, double *rb, double *x, double *y)
Definition: intersect.c:80
#define D1
Definition: intersect.c:75
#define SWAP(x, y)
Definition: intersect.c:78
#define D2
Definition: intersect.c:76
#define D
Definition: intersect.c:74
#define x