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