linecros.c

Go to the documentation of this file.

/*
 ****************************************************************************
 *
 * MODULE:       Vector library 
 *              
 * AUTHOR(S):    Original author CERL, probably Dave Gerdes.
 *               Update to GRASS 5.7 Radim Blazek.
 *
 * PURPOSE:      Lower level functions for reading/writing/manipulating vectors.
 *
 * COPYRIGHT:    (C) 2001 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *              License (>=v2). Read the file COPYING that comes with GRASS
 *              for details.
 *
 *****************************************************************************/
#include <stdio.h>

/***************************************************************
* test_for_intersection (ax1,ay1,ax2,ay2,bx1,by1,bx2,by2)
*   double ax1,ax2,ay1,ay2;
*   double bx1,bx2,by1,by2;
*
* returns
*   0 no intersection at all
*   1 the line segments intersect at only one point
*  -1 the line segments intersect at many points, i.e., overlapping
*     segments from the same line
*
* find_intersection (ax1,ay1,ax2,ay2,bx1,by1,bx2,by2,x,y)
*   double ax1,ax2,ay1,ay2;
*   double bx1,bx2,by1,by2;
*   double *x,*y;
*
* returns
*   0 no intersection
*   1 x,y set to (unique) intersection
*  -1 lines overlap, no unique intersection
*
* Based on the following:
*
*    (ax2-ax1)r1 - (bx2-bx1)r2 = ax2 - ax1
*    (ay2-ay1)r1 - (by2-by1)r2 = ay2 - ay1
*
* Solving for r1 and r2, if r1 and r2 are between 0 and 1,
* then line segments (ax1,ay1)(ax2,ay2) and (bx1,by1)(bx2,by2)
* intersect
****************************************************************/

#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))

int
dig_test_for_intersection(double ax1, double ay1,
                          double ax2, double ay2,
                          double bx1, double by1, double bx2, double by2)
{
    register double d, d1, d2;
    double t;
    int switched;

    if (ax1 > ax2 || (ax1 == ax2 && ay1 > ay2)) {
        t = ax1;
        ax1 = ax2;
        ax2 = t;

        t = ay1;
        ay1 = ay2;
        ay2 = t;
    }

    if (bx1 > bx2 || (bx1 == bx2 && by1 > by2)) {
        t = bx1;
        bx1 = bx2;
        bx2 = t;

        t = by1;
        by1 = by2;
        by2 = t;
    }

    switched = 0;
    if (bx1 < ax1)
        switched = 1;
    else if (bx1 == ax1) {
        if (bx2 < ax2)
            switched = 1;
        else if (bx2 == ax2) {
            if (by1 < ay1)
                switched = 1;
            else if (by1 == ay1) {
                if (by2 < ay2)
                    switched = 1;
            }
        }
    }

    if (switched) {
        t = ax1;
        ax1 = bx1;
        bx1 = t;
        t = ax2;
        ax2 = bx2;
        bx2 = t;

        t = ay1;
        ay1 = by1;
        by1 = t;
        t = ay2;
        ay2 = by2;
        by2 = t;
    }

    d = D;
    d1 = D1;
    d2 = D2;

    if (d > 0)
        return (d1 >= 0 && d2 >= 0 && d >= d1 && d >= d2);
    if (d < 0)
        return (d1 <= 0 && d2 <= 0 && d <= d1 && d <= d2);

    /* lines are parallel */
    if (d1 || d2)
        return 0;

    /* segments are colinear. check for overlap */

    /* Collinear vertical */
    if (ax1 == ax2) {
        if (ay1 > ay2) {
            t = ay1;
            ay1 = ay2;
            ay2 = t;
        }
        if (by1 > by2) {
            t = by1;
            by1 = by2;
            by2 = t;
        }
        if (ay1 > by2)
            return 0;
        if (ay2 < by1)
            return 0;

        /* there is overlap */

        if (ay1 == by2 || ay2 == by1)
            return 1;           /* endpoints only */

        return -1;                      /* true overlap   */
    }
    else {
        if (ax1 > ax2) {
            t = ax1;
            ax1 = ax2;
            ax2 = t;
        }
        if (bx1 > bx2) {
            t = bx1;
            bx1 = bx2;
            bx2 = t;
        }
        if (ax1 > bx2)
            return 0;
        if (ax2 < bx1)
            return 0;

        /* there is overlap */

        if (ax1 == bx2 || ax2 == bx1)
            return 1;           /* endpoints only */

        return -1;                      /* true overlap   */
    }    
    return 0; /* should not be reached */
}


int
dig_find_intersection(double ax1, double ay1,
                      double ax2, double ay2,
                      double bx1, double by1,
                      double bx2, double by2, double *x, double *y)
{
    register double d, r1, r2;
    double t;
    int switched;

    if (ax1 > ax2 || (ax1 == ax2 && ay1 > ay2)) {
        t = ax1;
        ax1 = ax2;
        ax2 = t;

        t = ay1;
        ay1 = ay2;
        ay2 = t;
    }

    if (bx1 > bx2 || (bx1 == bx2 && by1 > by2)) {
        t = bx1;
        bx1 = bx2;
        bx2 = t;

        t = by1;
        by1 = by2;
        by2 = t;
    }

    switched = 0;
    if (bx1 < ax1)
        switched = 1;
    else if (bx1 == ax1) {
        if (bx2 < ax2)
            switched = 1;
        else if (bx2 == ax2) {
            if (by1 < ay1)
                switched = 1;
            else if (by1 == ay1) {
                if (by2 < ay2)
                    switched = 1;
            }
        }
    }

    if (switched) {
        t = ax1;
        ax1 = bx1;
        bx1 = t;
        t = ax2;
        ax2 = bx2;
        bx2 = t;

        t = ay1;
        ay1 = by1;
        by1 = t;
        t = ay2;
        ay2 = by2;
        by2 = t;
    }

    d = D;

    if (d) {

        r1 = D1 / d;
        r2 = D2 / d;
        if (r1 < 0 || r1 > 1 || r2 < 0 || r2 > 1) {
            return 0;
        }
        *x = ax1 + r1 * (ax2 - ax1);
        *y = ay1 + r1 * (ay2 - ay1);
        return 1;
    }

    /* lines are parallel */
    if (D1 || D2) {
        return 0;
    }

    /* segments are colinear. check for overlap */

    /* Collinear vertical */
    if (ax1 == ax2) {
        if (ay1 > by2)
            return 0;
        if (ay2 < by1)
            return 0;

        /* 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 -1;
    }
    else {
        if (ax1 > bx2)
            return 0;
        if (ax2 < bx1)
            return 0;

        /* 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 -1;
    }

    return 0; /* should not be reached */
}