GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-e0f91531fb
gis/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 
84 int 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: gis/intersect.c:84
#define D1
Definition: gis/intersect.c:73
#define SWAP(x, y)
Definition: gis/intersect.c:76
#define D2
Definition: gis/intersect.c:74
#define D
Definition: gis/intersect.c:72
#define x