GRASS GIS 7 Programmer's Manual
7.8.4(2020)-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
80
int
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
}
G_intersect_line_segments
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
x
#define x
D1
#define D1
Definition:
intersect.c:75
SWAP
#define SWAP(x, y)
Definition:
intersect.c:78
D
#define D
Definition:
intersect.c:74
D2
#define D2
Definition:
intersect.c:76
gis
intersect.c
Generated on Mon Oct 5 2020 08:56:03 for GRASS GIS 7 Programmer's Manual by
1.8.18