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