GRASS GIS 7 Programmer's Manual  7.7.svn(2018)-r73666
gis/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 }
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:80
#define D2
Definition: gis/intersect.c:76
#define x
#define D1
Definition: gis/intersect.c:75
#define D
Definition: gis/intersect.c:74
#define SWAP(x, y)
Definition: gis/intersect.c:78