gis/intersect.c

Go to the documentation of this file.

/**************************************************************
 * find intersection between two lines defined by points on the lines
 * line segment A is (ax1,ay1) to (ax2,ay2)
 * line segment B is (bx1,by1) to (bx2,by2)
 * returns
 *   -1 segment A and B do not intersect (parallel without overlap)
 *    0 segment A and B do not intersect but extensions do intersect
 *    1 intersection is a single point
 *    2 intersection is a line segment (colinear with overlap)
 * x,y intersection point
 * ra - ratio that the intersection divides A
 * rb - ratio that the intersection divides B
 *
 *                              B2
 *                              /
 *                             /
 *   r=p/(p+q) : A1---p-------*--q------A2
 *                           /
 *                          /
 *                         B1
 *
 **************************************************************/
 
/**************************************************************
 *
 * A point P which lies on line defined by points A1=(x1,y1) and A2=(x2,y2)
 * is given by the equation r * (x2,y2) + (1-r) * (x1,y1).
 * if r is between 0 and 1, p lies between A1 and A2.
 *
 * Suppose points on line (A1, A2) has equation
 *     (x,y) = ra * (ax2,ay2) + (1-ra) * (ax1,ay1)
 * or for x and y separately
 *     x = ra * ax2 - ra * ax1 + ax1
 *     y = ra * ay2 - ra * ay1 + ay1
 * and the points on line (B1, B2) are represented by
 *     (x,y) = rb * (bx2,by2) + (1-rb) * (bx1,by1)
 * or for x and y separately
 *     x = rb * bx2 - rb * bx1 + bx1
 *     y = rb * by2 - rb * by1 + by1
 *
 * when the lines intersect, the point (x,y) has to
 * satisfy a system of 2 equations:
 *     ra * ax2 - ra * ax1 + ax1 = rb * bx2 - rb * bx1 + bx1
 *     ra * ay2 - ra * ay1 + ay1 = rb * by2 - rb * by1 + by1
 *
 * or
 *
 *     (ax2 - ax1) * ra - (bx2 - bx1) * rb = bx1 - ax1
 *     (ay2 - ay1) * ra - (by2 - by1) * rb = by1 - ay1
 *
 * by Cramer's method, one can solve this by computing 3
 * determinants of matrices:
 *
 *    M  = (ax2-ax1)  (bx1-bx2)
 *         (ay2-ay1)  (by1-by2)
 *
 *    M1 = (bx1-ax1)  (bx1-bx2)
 *         (by1-ay1)  (by1-by2)
 *
 *    M2 = (ax2-ax1)  (bx1-ax1)
 *         (ay2-ay1)  (by1-ay1)
 *
 * Which are exactly the determinants D, D2, D1 below:
 *
 *   D  ((ax2-ax1)*(by1-by2) - (ay2-ay1)*(bx1-bx2))
 *
 *   D1 ((bx1-ax1)*(by1-by2) - (by1-ay1)*(bx1-bx2))
 *
 *   D2 ((ax2-ax1)*(by1-ay1) - (ay2-ay1)*(bx1-ax1))
 ***********************************************************************/
 
#define D  ((ax2 - ax1) * (by1 - by2) - (ay2 - ay1) * (bx1 - bx2))
#define D1 ((bx1 - ax1) * (by1 - by2) - (by1 - ay1) * (bx1 - bx2))
#define D2 ((ax2 - ax1) * (by1 - ay1) - (ay2 - ay1) * (bx1 - ax1))
 
#define SWAP(x, y) \
    {              \
        double t;  \
        t = x;     \
        x = y;     \
        y = t;     \
    }
 
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)
{
    double d;
 
    if (ax1 > ax2 || (ax1 == ax2 && ay1 > ay2)) {
        SWAP(ax1, ax2)
        SWAP(ay1, ay2)
    }
 
    if (bx1 > bx2 || (bx1 == bx2 && by1 > by2)) {
        SWAP(bx1, bx2)
        SWAP(by1, by2)
    }
 
    d = D;
 
    if (d) { /* lines are not parallel */
        *ra = D1 / d;
        *rb = D2 / d;
 
        *x = ax1 + (*ra) * (ax2 - ax1);
        *y = ay1 + (*ra) * (ay2 - ay1);
        return (*ra >= 0.0 && *ra <= 1.0 && *rb >= 0.0 && *rb <= 1.0);
    }
 
    /* lines are parallel */
    if (D1 || D2) {
        return -1;
    }
 
    /* segments are colinear. check for overlap */
 
    /* Collinear vertical */
    if (ax1 == ax2) {
        if (ay1 > by2) {
            *x = ax1;
            *y = ay1;
            return 0; /* extensions overlap */
        }
        if (ay2 < by1) {
            *x = ax2;
            *y = ay2;
            return 0; /* extensions overlap */
        }
 
        /* there is overlap */
 
        if (ay1 == by2) {
            *x = ax1;
            *y = ay1;
            return 1; /* endpoints only */
        }
        if (ay2 == by1) {
            *x = ax2;
            *y = ay2;
            return 1; /* endpoints only */
        }
 
        /* overlap, no single intersection point */
        if (ay1 > by1 && ay1 < by2) {
            *x = ax1;
            *y = ay1;
        }
        else {
            *x = ax2;
            *y = ay2;
        }
        return 2;
    }
    else {
        if (ax1 > bx2) {
            *x = ax1;
            *y = ay1;
            return 0; /* extensions overlap */
        }
        if (ax2 < bx1) {
            *x = ax2;
            *y = ay2;
            return 0; /* extensions overlap */
        }
 
        /* there is overlap */
 
        if (ax1 == bx2) {
            *x = ax1;
            *y = ay1;
            return 1; /* endpoints only */
        }
        if (ax2 == bx1) {
            *x = ax2;
            *y = ay2;
            return 1; /* endpoints only */
        }
 
        /* overlap, no single intersection point */
        if (ax1 > bx1 && ax1 < bx2) {
            *x = ax1;
            *y = ay1;
        }
        else {
            *x = ax2;
            *y = ay2;
        }
        return 2;
    }
 
    return 0; /* should not be reached */
}