GRASS Programmer's Manual
6.4.1(2011)
|
00001 00002 /************************************************************** 00003 * find interesection between two lines defined by points on the lines 00004 * line segment A is (ax1,ay1) to (ax2,ay2) 00005 * line segment B is (bx1,by1) to (bx2,by2) 00006 * returns 00007 * -1 segment A and B do not intersect (parallel without overlap) 00008 * 0 segment A and B do not intersect but extensions do intersect 00009 * 1 intersection is a single point 00010 * 2 intersection is a line segment (colinear with overlap) 00011 * x,y intersection point 00012 * ra - ratio that the intersection divides A 00013 * rb - ratio that the intersection divides B 00014 * 00015 * B2 00016 * / 00017 * / 00018 * r=p/(p+q) : A1---p-------*--q------A2 00019 * / 00020 * / 00021 * B1 00022 * 00023 **************************************************************/ 00024 00025 /************************************************************** 00026 * 00027 * A point P which lies on line defined by points A1=(x1,y1) and A2=(x2,y2) 00028 * is given by the equation r * (x2,y2) + (1-r) * (x1,y1). 00029 * if r is between 0 and 1, p lies between A1 and A2. 00030 * 00031 * Suppose points on line (A1, A2) has equation 00032 * (x,y) = ra * (ax2,ay2) + (1-ra) * (ax1,ay1) 00033 * or for x and y separately 00034 * x = ra * ax2 - ra * ax1 + ax1 00035 * y = ra * ay2 - ra * ay1 + ay1 00036 * and the points on line (B1, B2) are represented by 00037 * (x,y) = rb * (bx2,by2) + (1-rb) * (bx1,by1) 00038 * or for x and y separately 00039 * x = rb * bx2 - rb * bx1 + bx1 00040 * y = rb * by2 - rb * by1 + by1 00041 * 00042 * when the lines intersect, the point (x,y) has to 00043 * satisfy a system of 2 equations: 00044 * ra * ax2 - ra * ax1 + ax1 = rb * bx2 - rb * bx1 + bx1 00045 * ra * ay2 - ra * ay1 + ay1 = rb * by2 - rb * by1 + by1 00046 * 00047 * or 00048 * 00049 * (ax2 - ax1) * ra - (bx2 - bx1) * rb = bx1 - ax1 00050 * (ay2 - ay1) * ra - (by2 - by1) * rb = by1 - ay1 00051 * 00052 * by Cramer's method, one can solve this by computing 3 00053 * determinants of matrices: 00054 * 00055 * M = (ax2-ax1) (bx1-bx2) 00056 * (ay2-ay1) (by1-by2) 00057 * 00058 * M1 = (bx1-ax1) (bx1-bx2) 00059 * (by1-ay1) (by1-by2) 00060 * 00061 * M2 = (ax2-ax1) (bx1-ax1) 00062 * (ay2-ay1) (by1-ay1) 00063 * 00064 * Which are exactly the determinants D, D2, D1 below: 00065 * 00066 * D ((ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2)) 00067 * 00068 * D1 ((bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2)) 00069 * 00070 * D2 ((ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1)) 00071 ***********************************************************************/ 00072 00073 00074 #define D ((ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2)) 00075 #define D1 ((bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2)) 00076 #define D2 ((ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1)) 00077 00078 #define SWAP(x,y) {int t; t=x; x=y; y=t;} 00079 00080 int G_intersect_line_segments(double ax1, double ay1, double ax2, double ay2, 00081 double bx1, double by1, double bx2, double by2, 00082 double *ra, double *rb, double *x, double *y) 00083 { 00084 double d; 00085 00086 d = D; 00087 00088 if (d) { /* lines are not parallel */ 00089 *ra = D1 / d; 00090 *rb = D2 / d; 00091 00092 *x = ax1 + (*ra) * (ax2 - ax1); 00093 *y = ay1 + (*ra) * (ay2 - ay1); 00094 return (*ra >= 0.0 && *ra <= 1.0 && *rb >= 0.0 && *rb <= 1.0); 00095 } 00096 00097 if (D1 || D2) 00098 return -1; /* lines are parallel, not colinear */ 00099 00100 if (ax1 > ax2) { 00101 SWAP(ax1, ax2) 00102 } 00103 if (bx1 > bx2) { 00104 SWAP(bx1, bx2) 00105 } 00106 if (ax1 > bx2) 00107 return -1; 00108 if (ax2 < bx1) 00109 return -1; 00110 00111 /* there is overlap */ 00112 if (ax1 == bx2) { 00113 *x = ax1; 00114 *y = ay1; 00115 return 1; /* at endpoints only */ 00116 } 00117 if (ax2 == bx1) { 00118 *x = ax2; 00119 *y = ay2; 00120 return 1; /* at endpoints only */ 00121 } 00122 00123 return 2; /* colinear with overlap on an interval, not just a single point */ 00124 }