intersect.c

Go to the documentation of this file.
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 (
00081     double ax1,double ay1, double ax2,double ay2,
00082     double bx1,double by1, double bx2,double by2,
00083     double *ra,double *rb,
00084     double *x,double *y)
00085 {
00086     double d;
00087 
00088     d = D;
00089 
00090     if (d) /* lines are not parallel */
00091     {
00092         *ra = D1/d;
00093         *rb = D2/d;
00094 
00095         *x = ax1 + (*ra) * (ax2 - ax1) ;
00096         *y = ay1 + (*ra) * (ay2 - ay1) ;
00097         return (*ra >= 0.0 && *ra <= 1.0 && *rb >= 0.0 && *rb <= 1.0);
00098     }
00099 
00100     if (D1 || D2) return -1;  /* lines are parallel, not colinear */
00101 
00102     if (ax1 > ax2)
00103     {
00104         SWAP (ax1, ax2)
00105     }
00106     if (bx1 > bx2)
00107     {
00108         SWAP (bx1, bx2)
00109     }
00110     if (ax1 > bx2) return -1;
00111     if (ax2 < bx1) return -1;
00112 
00113     /* there is overlap */
00114     if (ax1 == bx2)
00115     {
00116         *x = ax1;
00117         *y = ay1;
00118         return 1; /* at endpoints only */
00119     }
00120     if (ax2 == bx1)
00121     {
00122         *x = ax2;
00123         *y = ay2;
00124         return 1; /* at endpoints only */
00125     }
00126 
00127     return 2; /* colinear with overlap on an interval, not just a single point*/
00128 }

Generated on Wed Dec 19 14:59:06 2007 for GRASS by  doxygen 1.5.4