pointpoly.c

Go to the documentation of this file.

/*
Here it is, code for the angle test, the barycentric test, and the two
crossings tests, with a main program to test their speeds.  You may need to
use something other than srand()/drand48() for your random number generator,
and the times() command for timing.  You'll also want to change the
*_TEST_TIMES constants if you're using something slower than an HP 720
workstation.  Feel free to complain about the slowness of any of the code -
I'm always interested in new hacks.
 
 
 * Point in polygon inside/outside tester.  Angle summation, barycentric
 * coordinates, and ray along x-axis (crossings testing) compared.
 *
 * copyright 1992 by Eric Haines, 3D/Eye Inc, erich@eye.com
 */
 
#include <math.h>
#include <sys/types.h>
#include <sys/param.h>
#include <sys/times.h>
 
#define X   0
#define Y   1
 
double  drand48() ;
 
double  AngleTimeTotal ;
double  BarycentricTimeTotal ;
double  CrossingsTimeTotal ;
double  MacmartinTimeTotal ;
 
/* minimum & maximum number of polygon vertices to generate */
#define MIN_VERTS   3
#define MAX_VERTS   7
 
/* number of different polygons to try */
#define TEST_POLYGONS   30
 
/* number of different intersection points to try */
#define TEST_POINTS 30
 
/* number of times to try a single point vs. a polygon */
/* this should be > 1 / ( HZ * approx. single test time in seconds ) */
#define ANGLE_TEST_TIMES    1000
#define BARYCENTRIC_TEST_TIMES  10000
#define CROSSINGS_TEST_TIMES    10000
#define MACMARTIN_TEST_TIMES    10000
 
main(argc,argv)
int argc;  char *argv[];
{
int i, j, numverts, inside_tot ;
int angle_flag, barycentric_flag, crossings_flag, macmartin_flag ;
double  pgon[MAX_VERTS][2] ;
double  point[2] ;
 
    srand( 12345 ) ;
 
    AngleTimeTotal = 0.0 ;
    BarycentricTimeTotal = 0.0 ;
    CrossingsTimeTotal = 0.0 ;
    MacmartinTimeTotal = 0.0 ;
    inside_tot = 0 ;
 
    for ( i = 0 ; i < TEST_POLYGONS ; i++ ) {
 
#ifdef CENTERED_SQUARE
    /* for debugging purposes, test against a square */
    numverts = 4 ;
    pgon[0][X] = 0.2 ;
    pgon[0][Y] = 0.2 ;
    pgon[1][X] = 0.8 ;
    pgon[1][Y] = 0.2 ;
    pgon[2][X] = 0.8 ;
    pgon[2][Y] = 0.8 ;
    pgon[3][X] = 0.2 ;
    pgon[3][Y] = 0.8 ;
#else
    /* make an arbitrary polygon fitting 0-1 range in x and y */
    numverts = MIN_VERTS +
        (int)(drand48() * (double)(MAX_VERTS-MIN_VERTS+1)) ;
    for ( j = 0 ; j < numverts ; j++ ) {
        pgon[j][X] = drand48() ;
        pgon[j][Y] = drand48() ;
    }
#endif
 
    /* now try # of points against it */
    for ( j = 0 ; j < TEST_POINTS ; j++ ) {
        point[X] = drand48() ;
        point[Y] = drand48() ;
        angle_flag = angletest( pgon, numverts, point ) ;
        barycentric_flag = barycentrictest( pgon, numverts, point ) ;
        crossings_flag = crossingstest( pgon, numverts, point ) ;
        macmartin_flag = macmartintest( pgon, numverts, point ) ;
 
        /* reality check */
        if ( angle_flag != crossings_flag ) {
        printf( "angle test says %s, crossings test says %s\n",
            angle_flag ? "INSIDE" : "OUTSIDE",
            crossings_flag ? "INSIDE" : "OUTSIDE" ) ;
        printf( "point %g %g\n", (float)point[X], (float)point[Y] ) ;
        printf( "polygon:\n" ) ;
        for ( j = 0 ; j < numverts ; j++ ) {
            printf( "  %g %g\n", (float)pgon[j][X], (float)pgon[j][Y]);
        }
        }
        if ( barycentric_flag != macmartin_flag ) {
        printf( "barycentric test says %s, macmartin test says %s\n",
            barycentric_flag ? "INSIDE" : "OUTSIDE",
            macmartin_flag ? "INSIDE" : "OUTSIDE" ) ;
        printf( "point %g %g\n", (float)point[X], (float)point[Y] ) ;
        printf( "polygon:\n" ) ;
        for ( j = 0 ; j < numverts ; j++ ) {
            printf( "  %g %g\n", (float)pgon[j][X], (float)pgon[j][Y]);
        }
        }
 
        inside_tot += crossings_flag ;
    }
    }
    printf( "angle test time: %g microseconds per test\n",
    (float)( AngleTimeTotal/(double)(TEST_POINTS*TEST_POLYGONS) ) ) ;
    printf( "barycentric test time: %g microseconds per test\n",
    (float)( BarycentricTimeTotal/(double)(TEST_POINTS*TEST_POLYGONS) ) ) ;
    printf( "crossings test time: %g microseconds per test\n",
    (float)( CrossingsTimeTotal/(double)(TEST_POINTS*TEST_POLYGONS) ) ) ;
    printf( "macmartin crossings test time: %g microseconds per test\n",
    (float)( MacmartinTimeTotal/(double)(TEST_POINTS*TEST_POLYGONS) ) ) ;
 
    printf( "%g %% of all points were inside polygons\n",
    (float)inside_tot * 100.0 / (float)(TEST_POINTS*TEST_POLYGONS) ) ;
}
 
/* sum angles of vtxN - point - vtxN+1, check if in PI to 3*PI range */
int
angletest( pgon, numverts, point )
double  pgon[MAX_VERTS][2] ;
int numverts ;
double  point[2] ;
{
int i, j, inside_flag ;
struct  tms     timebuf ;
long    timestart ;
long    timestop ;
double  *vtx0, *vtx1, angle, len, vec0[2], vec1[2] ;
 
    /* do the test a bunch of times to get a useful time reading */
    timestart = times( &timebuf ) ;
    for ( i = 0 ; i < ANGLE_TEST_TIMES ; i++ ) {
    /* sum the angles and see if answer mod 2*PI > PI */
    vtx0 = pgon[numverts-1] ;
    vec0[X] = vtx0[X] - point[X] ;
    vec0[Y] = vtx0[Y] - point[Y] ;
    if ( (len = hypot( vec0[X], vec0[Y] )) <= 0.0 ) {
        /* point and vertex coincide */
        return( 1 ) ;
    }
    vec0[X] /= len ;
    vec0[Y] /= len ;
 
    angle = 0.0 ;
    for ( j = 0 ; j < numverts ; j++ ) {
        vtx1 = pgon[j] ;
        vec1[X] = vtx1[X] - point[X] ;
        vec1[Y] = vtx1[Y] - point[Y] ;
        if ( (len = hypot( vec1[X], vec1[Y] )) <= 0.0 ) {
        /* point and vertex coincide */
        return( 1 ) ;
        }
        vec1[X] /= len ;
        vec1[Y] /= len ;
 
        /* check if vec1 is to "left" or "right" of vec0 */
        if ( vec0[X] * vec1[Y] - vec1[X] * vec0[Y] >= 0.0 ) {
        /* add angle due to dot product of vectors */
        angle += acos( vec0[X] * vec1[X] + vec0[Y] * vec1[Y] ) ;
        } else {
        /* subtract angle due to dot product of vectors */
        angle -= acos( vec0[X] * vec1[X] + vec0[Y] * vec1[Y] ) ;
        }
 
        /* get to next point */
        vtx0 = vtx1 ;
        vec0[X] = vec1[X] ;
        vec0[Y] = vec1[Y] ;
    }
    /* test if between PI and 3*PI, 5*PI and 7*PI, etc */
    /* if we care about is winding number > 0, then just:
           inside_flag = fabs(angle) > M_PI ;
     */
    inside_flag = fmod( fabs(angle) + M_PI, 4.0*M_PI ) > 2.0*M_PI ;
    }
    timestop = times( &timebuf ) ;
    /* time in microseconds */
    AngleTimeTotal += 1000000.0 * (double)(timestop - timestart) /
    (double)(HZ * ANGLE_TEST_TIMES) ;
 
    return (inside_flag) ;
}
 
/* barycentric, a la Gems I, with a little efficiency tuning */
int
barycentrictest( pgon, numverts, point )
double  pgon[MAX_VERTS][2] ;
int numverts ;
double  point[2] ;
{
int i, tris_hit, inside_flag, p1, p2 ;
struct  tms     timebuf ;
long    timestart ;
long    timestop ;
double  tx, ty, u0, u1, u2, v0, v1, alpha, beta, denom ;
 
    /* do the test a bunch of times to get a useful time reading */
    timestart = times( &timebuf ) ;
    for ( i = 0 ; i < BARYCENTRIC_TEST_TIMES ; i++ ) {
 
    tx = point[X] ;
    ty = point[Y] ;
 
    tris_hit = 0 ;
 
    for ( p1 = 1, p2 = 2 ; p2 < numverts ; p1++, p2++ ) {
 
        if ( ( u1 = pgon[0][X] - pgon[p2][X] ) == 0.0 ) {
 
        /* zero area test - optional */
        if ( ( u2 = pgon[p1][X] - pgon[p2][X] ) == 0.0 ) {
            goto NextTri;
        }
 
        /* Compute intersection point */
        if ( ( ( beta = ( tx - pgon[p2][X] ) / u2 ) < 0.0 ) ||
               ( beta > 1.0 ) ) {
 
            goto NextTri;
        }
 
        if ( ( v1 = pgon[0][Y] - pgon[p2][Y] ) == 0.0 ) {
            goto NextTri;
        }
 
        if ( ( alpha = ( ty - pgon[p2][Y] - beta *
            ( pgon[p1][Y] - pgon[p2][Y] ) ) / v1 ) < 0.0 ) {
            goto NextTri;
        }
 
        } else {
        if ( ( denom = ( pgon[p1][Y] - pgon[p2][Y] ) * u1 -
            ( u2 = pgon[p1][X] - pgon[p2][X] ) *
            ( v1 = pgon[0][Y] - pgon[p2][Y] ) ) == 0.0 ) {
 
            goto NextTri;
        }
 
        /* Compute intersection point & subtract 0 vertex */
        u0 = tx - pgon[p2][X] ;
        v0 = ty - pgon[p2][Y] ;
 
        if ( ( ( ( beta = ( v0 * u1 - u0 * v1 ) / denom ) ) < 0.0 ) ||
               ( beta > 1.0 ) ) {
 
            goto NextTri;
        }
        if ( ( alpha = ( u0 - beta * u2 ) / u1 ) < 0.0 ) {
            goto NextTri;
        }
        }
 
        /* check gamma */
        if ( alpha + beta <= 1.0 ) {
        /* survived */
        tris_hit++ ;
        }
 
        NextTri: ;
    }
    inside_flag = tris_hit & 0x1 ;
    }
    timestop = times( &timebuf ) ;
    /* time in microseconds */
    BarycentricTimeTotal += 1000000.0 * (double)(timestop - timestart) /
    (double)(HZ * BARYCENTRIC_TEST_TIMES) ;
 
    return (inside_flag) ;
}
 
/* shoot a test ray along +X axis - slower version, less messy */
int
crossingstest( pgon, numverts, point )
double  pgon[MAX_VERTS][2] ;
int numverts ;
double  point[2] ;
{
int i, j, inside_flag, xflag0 ;
struct  tms     timebuf ;
long    timestart ;
long    timestop ;
double  *vtx0, *vtx1, dv0 ;
int crossings, yflag0, yflag1 ;
 
    /* do the test a bunch of times to get a useful time reading */
    timestart = times( &timebuf ) ;
    for ( i = 0 ; i < CROSSINGS_TEST_TIMES ; i++ ) {
 
    vtx0 = pgon[numverts-1] ;
    /* get test bit for above/below Y axis */
    yflag0 = ( dv0 = vtx0[Y] - point[Y] ) >= 0.0 ;
 
    crossings = 0 ;
    for ( j = 0 ; j < numverts ; j++ ) {
        /* cleverness:  bobble between filling endpoints of edges, so
         * that the previous edge's shared endpoint is maintained.
         */
        if ( j & 0x1 ) {
        vtx0 = pgon[j] ;
        yflag0 = ( dv0 = vtx0[Y] - point[Y] ) >= 0.0 ;
        } else {
        vtx1 = pgon[j] ;
        yflag1 = ( vtx1[Y] >= point[Y] ) ;
        }
 
        /* check if points not both above/below X axis - can't hit ray */
        if ( yflag0 != yflag1 ) {
        /* check if points on same side of Y axis */
        if ( ( xflag0 = ( vtx0[X] >= point[X] ) ) ==
             ( vtx1[X] >= point[X] ) ) {
 
            if ( xflag0 ) crossings++ ;
        } else {
            /* compute intersection of pgon segment with X ray, note
             * if > point's X.
             */
            crossings += (vtx0[X] -
            dv0*( vtx1[X]-vtx0[X])/(vtx1[Y]-vtx0[Y])) >= point[X] ;
        }
        }
    }
    /* test if crossings is odd */
    /* if we care about is winding number > 0, then just:
           inside_flag = crossings > 0 ;
     */
    inside_flag = crossings & 0x01 ;
    }
    timestop = times( &timebuf ) ;
    /* time in microseconds */
    CrossingsTimeTotal += 1000000.0 * (double)(timestop - timestart) /
    (double)(HZ * CROSSINGS_TEST_TIMES) ;
 
    return (inside_flag) ;
}
 
/* shoot a test ray along +X axis - macmartinized by me, a bit messier */
int
macmartintest( pgon, numverts, point )
double  pgon[MAX_VERTS][2] ;
int numverts ;
double  point[2] ;
{
int i, inside_flag, xflag0 ;
struct  tms     timebuf ;
long    timestart ;
long    timestop ;
double  *p, *stop ;
int crossings ;
double  tx, ty, y ;
 
    /* do the test a bunch of times to get a useful time reading */
    timestart = times( &timebuf ) ;
    for ( i = 0 ; i < MACMARTIN_TEST_TIMES ; i++ ) {
 
    crossings = 0 ;
 
    tx = point[X] ;
    ty = point[Y] ;
    y = pgon[numverts-1][Y] ;
 
    p = (double *)pgon + 1 ;
    if ( ( y >= ty ) != ( *p >= ty ) ) {
        /* x crossing */
        if ( ( xflag0 = ( pgon[numverts-1][X] >= tx ) ) ==
         ( *(double *)pgon >= tx ) ) {
 
        if ( xflag0 ) crossings++ ;
        } else {
        /* compute intersection of pgon segment with X ray, note
         * if > point's X.
         */
        crossings += ( pgon[numverts-1][X] -
        (y-ty)*( *(double *)pgon - pgon[numverts-1][X])/(*p-y)) >= tx ;
        }
    }
 
    stop = pgon[numverts] ;
 
    for ( y=*p,p += 2 ; p < stop ; y=*p,p+=2) {
 
        if ( y >= ty ) {
        while ( (p < stop) && (*p >= ty) ) p+=2 ;
        if ( p >= stop ) goto Exit ;
        /* y crosses */
        if ( ( xflag0 = ( *(p-3) >= tx ) ) ==
             ( *(p-1) >= tx ) ) {
 
            if ( xflag0 ) crossings++ ;
        } else {
            /* compute intersection of pgon segment with X ray, note
             * if > point's X.
             */
            crossings += ( *(p-3) -
            (*(p-2)-ty)*( *(p-1)-*(p-3))/(*p-*(p-2))) >= tx ;
        }
        } else {
        while ( (p < stop) && (*p < ty)) p+=2 ;
        if ( p >= stop ) goto Exit ;
        /* y crosses */
        if ( ( xflag0 = ( *(p-3) >= tx ) ) ==
             ( *(p-1) >= tx ) ) {
 
            if ( xflag0 ) crossings++ ;
        } else {
            /* compute intersection of pgon segment with X ray, note
             * if > point's X.
             */
            crossings += ( *(p-3) -
            (*(p-2)-ty)*( *(p-1)-*(p-3))/(*p-*(p-2))) >= tx ;
        }
        }
    }
 
    Exit:
    /* test if crossings is odd */
    /* if we care about is winding number > 0, then just:
           inside_flag = crossings > 0 ;
     */
    inside_flag = crossings & 0x01 ;
    }
    timestop = times( &timebuf ) ;
    /* time in microseconds */
    MacmartinTimeTotal += 1000000.0 * (double)(timestop - timestart) /
    (double)(HZ * MACMARTIN_TEST_TIMES) ;
 
    return (inside_flag) ;
}