GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-535c39c9fc
linecros.c
Go to the documentation of this file.
1 /*****************************************************************************
2  *
3  * MODULE: Vector library
4  *
5  * AUTHOR(S): Original author CERL, probably Dave Gerdes.
6  * Update to GRASS 5.7 Radim Blazek.
7  *
8  * PURPOSE: Lower level functions for reading/writing/manipulating vectors.
9  *
10  * COPYRIGHT: (C) 2001 by the GRASS Development Team
11  *
12  * This program is free software under the GNU General Public
13  * License (>=v2). Read the file COPYING that comes with GRASS
14  * for details.
15  *
16  *****************************************************************************/
17 
18 #include <stdio.h>
19 
20 /***************************************************************
21  * test_for_intersection (ax1,ay1,ax2,ay2,bx1,by1,bx2,by2)
22  * double ax1,ax2,ay1,ay2;
23  * double bx1,bx2,by1,by2;
24  *
25  * returns
26  * 0 no intersection at all
27  * 1 the line segments intersect at only one point
28  * -1 the line segments intersect at many points, i.e., overlapping
29  * segments from the same line
30  *
31  * find_intersection (ax1,ay1,ax2,ay2,bx1,by1,bx2,by2,x,y)
32  * double ax1,ax2,ay1,ay2;
33  * double bx1,bx2,by1,by2;
34  * double *x,*y;
35  *
36  * returns
37  * 0 no intersection
38  * 1 x,y set to (unique) intersection
39  * -1 lines overlap, no unique intersection
40  *
41  * Based on the following:
42  *
43  * (ax2-ax1)r1 - (bx2-bx1)r2 = ax2 - ax1
44  * (ay2-ay1)r1 - (by2-by1)r2 = ay2 - ay1
45  *
46  * Solving for r1 and r2, if r1 and r2 are between 0 and 1,
47  * then line segments (ax1,ay1)(ax2,ay2) and (bx1,by1)(bx2,by2)
48  * intersect
49  ****************************************************************/
50 
51 #define D ((ax2 - ax1) * (by1 - by2) - (ay2 - ay1) * (bx1 - bx2))
52 
53 #define D1 ((bx1 - ax1) * (by1 - by2) - (by1 - ay1) * (bx1 - bx2))
54 
55 #define D2 ((ax2 - ax1) * (by1 - ay1) - (ay2 - ay1) * (bx1 - ax1))
56 
57 int dig_test_for_intersection(double ax1, double ay1, double ax2, double ay2,
58  double bx1, double by1, double bx2, double by2)
59 {
60  register double d, d1, d2;
61  double t;
62  int switched;
63 
64  if (ax1 > ax2 || (ax1 == ax2 && ay1 > ay2)) {
65  t = ax1;
66  ax1 = ax2;
67  ax2 = t;
68 
69  t = ay1;
70  ay1 = ay2;
71  ay2 = t;
72  }
73 
74  if (bx1 > bx2 || (bx1 == bx2 && by1 > by2)) {
75  t = bx1;
76  bx1 = bx2;
77  bx2 = t;
78 
79  t = by1;
80  by1 = by2;
81  by2 = t;
82  }
83 
84  switched = 0;
85  if (bx1 < ax1)
86  switched = 1;
87  else if (bx1 == ax1) {
88  if (bx2 < ax2)
89  switched = 1;
90  else if (bx2 == ax2) {
91  if (by1 < ay1)
92  switched = 1;
93  else if (by1 == ay1) {
94  if (by2 < ay2)
95  switched = 1;
96  }
97  }
98  }
99 
100  if (switched) {
101  t = ax1;
102  ax1 = bx1;
103  bx1 = t;
104  t = ax2;
105  ax2 = bx2;
106  bx2 = t;
107 
108  t = ay1;
109  ay1 = by1;
110  by1 = t;
111  t = ay2;
112  ay2 = by2;
113  by2 = t;
114  }
115 
116  d = D;
117  d1 = D1;
118  d2 = D2;
119 
120  if (d > 0)
121  return (d1 >= 0 && d2 >= 0 && d >= d1 && d >= d2);
122  if (d < 0)
123  return (d1 <= 0 && d2 <= 0 && d <= d1 && d <= d2);
124 
125  /* lines are parallel */
126  if (d1 || d2)
127  return 0;
128 
129  /* segments are colinear. check for overlap */
130 
131  /* Collinear vertical */
132  if (ax1 == ax2) {
133  if (ay1 > ay2) {
134  t = ay1;
135  ay1 = ay2;
136  ay2 = t;
137  }
138  if (by1 > by2) {
139  t = by1;
140  by1 = by2;
141  by2 = t;
142  }
143  if (ay1 > by2)
144  return 0;
145  if (ay2 < by1)
146  return 0;
147 
148  /* there is overlap */
149 
150  if (ay1 == by2 || ay2 == by1)
151  return 1; /* endpoints only */
152 
153  return -1; /* true overlap */
154  }
155  else {
156  if (ax1 > ax2) {
157  t = ax1;
158  ax1 = ax2;
159  ax2 = t;
160  }
161  if (bx1 > bx2) {
162  t = bx1;
163  bx1 = bx2;
164  bx2 = t;
165  }
166  if (ax1 > bx2)
167  return 0;
168  if (ax2 < bx1)
169  return 0;
170 
171  /* there is overlap */
172 
173  if (ax1 == bx2 || ax2 == bx1)
174  return 1; /* endpoints only */
175 
176  return -1; /* true overlap */
177  }
178  return 0; /* should not be reached */
179 }
180 
181 int dig_find_intersection(double ax1, double ay1, double ax2, double ay2,
182  double bx1, double by1, double bx2, double by2,
183  double *x, double *y)
184 {
185  register double d, r1, r2;
186  double t;
187  int switched;
188 
189  if (ax1 > ax2 || (ax1 == ax2 && ay1 > ay2)) {
190  t = ax1;
191  ax1 = ax2;
192  ax2 = t;
193 
194  t = ay1;
195  ay1 = ay2;
196  ay2 = t;
197  }
198 
199  if (bx1 > bx2 || (bx1 == bx2 && by1 > by2)) {
200  t = bx1;
201  bx1 = bx2;
202  bx2 = t;
203 
204  t = by1;
205  by1 = by2;
206  by2 = t;
207  }
208 
209  switched = 0;
210  if (bx1 < ax1)
211  switched = 1;
212  else if (bx1 == ax1) {
213  if (bx2 < ax2)
214  switched = 1;
215  else if (bx2 == ax2) {
216  if (by1 < ay1)
217  switched = 1;
218  else if (by1 == ay1) {
219  if (by2 < ay2)
220  switched = 1;
221  }
222  }
223  }
224 
225  if (switched) {
226  t = ax1;
227  ax1 = bx1;
228  bx1 = t;
229  t = ax2;
230  ax2 = bx2;
231  bx2 = t;
232 
233  t = ay1;
234  ay1 = by1;
235  by1 = t;
236  t = ay2;
237  ay2 = by2;
238  by2 = t;
239  }
240 
241  d = D;
242 
243  if (d) {
244 
245  r1 = D1 / d;
246  r2 = D2 / d;
247  if (r1 < 0 || r1 > 1 || r2 < 0 || r2 > 1) {
248  return 0;
249  }
250  *x = ax1 + r1 * (ax2 - ax1);
251  *y = ay1 + r1 * (ay2 - ay1);
252  return 1;
253  }
254 
255  /* lines are parallel */
256  if (D1 || D2) {
257  return 0;
258  }
259 
260  /* segments are colinear. check for overlap */
261 
262  /* Collinear vertical */
263  if (ax1 == ax2) {
264  if (ay1 > by2)
265  return 0;
266  if (ay2 < by1)
267  return 0;
268 
269  /* there is overlap */
270 
271  if (ay1 == by2) {
272  *x = ax1;
273  *y = ay1;
274  return 1; /* endpoints only */
275  }
276  if (ay2 == by1) {
277  *x = ax2;
278  *y = ay2;
279  return 1; /* endpoints only */
280  }
281 
282  /* overlap, no single intersection point */
283  if (ay1 > by1 && ay1 < by2) {
284  *x = ax1;
285  *y = ay1;
286  }
287  else {
288  *x = ax2;
289  *y = ay2;
290  }
291  return -1;
292  }
293  else {
294  if (ax1 > bx2)
295  return 0;
296  if (ax2 < bx1)
297  return 0;
298 
299  /* there is overlap */
300 
301  if (ax1 == bx2) {
302  *x = ax1;
303  *y = ay1;
304  return 1; /* endpoints only */
305  }
306  if (ax2 == bx1) {
307  *x = ax2;
308  *y = ay2;
309  return 1; /* endpoints only */
310  }
311 
312  /* overlap, no single intersection point */
313  if (ax1 > bx1 && ax1 < bx2) {
314  *x = ax1;
315  *y = ay1;
316  }
317  else {
318  *x = ax2;
319  *y = ay2;
320  }
321  return -1;
322  }
323 
324  return 0; /* should not be reached */
325 }
int dig_find_intersection(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2, double *x, double *y)
Definition: linecros.c:181
#define D1
Definition: linecros.c:53
#define D2
Definition: linecros.c:55
#define D
Definition: linecros.c:51
int dig_test_for_intersection(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2)
Definition: linecros.c:57
double t
Definition: r_raster.c:39
#define x