ptinpoly.c

[Game.Programming].Academic - Graphics Gems (6 books source code)

/* ptinpoly.c - point in polygon inside/outside code.

   by Eric Haines, 3D/Eye Inc, [e-mail]

   This code contains the following algorithms:
	crossings - count the crossing made by a ray from the test point
	crossings-multiply - as above, but avoids a division; often a bit faster
	angle summation - sum the angle formed by point and vertex pairs
	weiler angle summation - sum the angles using quad movements
	half-plane testing - test triangle fan using half-space planes
	barycentric coordinates - test triangle fan w/barycentric coords
	spackman barycentric - preprocessed barycentric coordinates
	trapezoid testing - bin sorting algorithm
	grid testing - grid imposed on polygon
	exterior test - for convex polygons, check exterior of polygon
	inclusion test - for convex polygons, use binary search for edge.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "ptinpoly.h"

#define X	0
#define Y	1

#ifndef TRUE
#define TRUE	1
#define FALSE	0
#endif

#ifndef HUGE
#define HUGE	1.797693134862315e+308
#endif

#ifndef M_PI
#define M_PI	3.14159265358979323846
#endif

/* test if a & b are within epsilon.  Favors cases where a < b */
#define Near(a,b,eps)	( ((b)-(eps)<(a)) && ((a)-(eps)<(b)) )

#define MALLOC_CHECK( a )	if ( !(a) ) {				   \
				    fprintf( stderr, "out of memory\n" ) ; \
				    exit(1) ;				   \
				}


/* ======= Crossings algorithm ============================================ */

/* Shoot a test ray along +X axis.  The strategy, from MacMartin, is to
 * compare vertex Y values to the testing point's Y and quickly discard
 * edges which are entirely to one side of the test ray.
 *
 * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
 * _point_, returns 1 if inside, 0 if outside.	WINDING and CONVEX can be
 * defined for this test.
 */
int CrossingsTest( pgon, numverts, point )
double	pgon[][2] ;
int	numverts ;
double	point[2] ;
{
#ifdef	WINDING
register int	crossings ;
#endif
register int	j, yflag0, yflag1, inside_flag, xflag0 ;
register double ty, tx, *vtx0, *vtx1 ;
#ifdef	CONVEX
register int	line_flag ;
#endif

    tx = point[X] ;
    ty = point[Y] ;

    vtx0 = pgon[numverts-1] ;
    /* get test bit for above/below X axis */
    yflag0 = ( vtx0[Y] >= ty ) ;
    vtx1 = pgon[0] ;

#ifdef	WINDING
    crossings = 0 ;
#else
    inside_flag = 0 ;
#endif
#ifdef	CONVEX
    line_flag = 0 ;
#endif
    for ( j = numverts+1 ; --j ; ) {

	yflag1 = ( vtx1[Y] >= ty ) ;
	/* check if endpoints straddle (are on opposite sides) of X axis
	 * (i.e. the Y's differ); if so, +X ray could intersect this edge.
	 */
	if ( yflag0 != yflag1 ) {
	    xflag0 = ( vtx0[X] >= tx ) ;
	    /* check if endpoints are on same side of the Y axis (i.e. X's
	     * are the same); if so, it's easy to test if edge hits or misses.
	     */
	    if ( xflag0 == ( vtx1[X] >= tx ) ) {

		/* if edge's X values both right of the point, must hit */
#ifdef	WINDING
		if ( xflag0 ) crossings += ( yflag0 ? -1 : 1 ) ;
#else
		if ( xflag0 ) inside_flag = !inside_flag ;
#endif
	    } else {
		/* compute intersection of pgon segment with +X ray, note
		 * if >= point's X; if so, the ray hits it.
		 */
		if ( (vtx1[X] - (vtx1[Y]-ty)*
		     ( vtx0[X]-vtx1[X])/(vtx0[Y]-vtx1[Y])) >= tx ) {
#ifdef	WINDING
		    crossings += ( yflag0 ? -1 : 1 ) ;
#else
		    inside_flag = !inside_flag ;
#endif
		}
	    }
#ifdef	CONVEX
	    /* if this is second edge hit, then done testing */
	    if ( line_flag ) goto Exit ;

	    /* note that one edge has been hit by the ray's line */
	    line_flag = TRUE ;
#endif
	}

	/* move to next pair of vertices, retaining info as possible */
	yflag0 = yflag1 ;
	vtx0 = vtx1 ;
	vtx1 += 2 ;
    }
#ifdef	CONVEX
    Exit: ;
#endif
#ifdef	WINDING
    /* test if crossings is not zero */
    inside_flag = (crossings != 0) ;
#endif

    return( inside_flag ) ;
}

/* ======= Angle summation algorithm ======================================= */

/* Sum angles made by (vtxN to test point to vtxN+1), check if in proper
 * range to be inside.	VERY SLOW, included for tutorial reasons (i.e.
 * to show why one should never use this algorithm).
 *
 * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
 * _point_, returns 1 if inside, 0 if outside.
 */
int AngleTest( pgon, numverts, point )
double	pgon[][2] ;
int	numverts ;
double	point[2] ;
{
register double *vtx0, *vtx1, angle, len, vec0[2], vec1[2], vec_dot ;
register int	j ;
int	inside_flag ;

    /* 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] ;
    len = sqrt( vec0[X] * vec0[X] + vec0[Y] * vec0[Y] ) ;
    if ( len <= 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] ;
	len = sqrt( vec1[X] * vec1[X] + vec1[Y] * vec1[Y] ) ;
	if ( len <= 0.0 ) {
	    /* point and vertex coincide */
	    return( 1 ) ;
	}
	vec1[X] /= len ;
	vec1[Y] /= len ;

	/* check if vec1 is to "left" or "right" of vec0 */
	vec_dot = vec0[X] * vec1[X] + vec0[Y] * vec1[Y] ;
	if ( vec_dot < -1.0 ) {
	    /* point is on edge, so always add 180 degrees.  always
	     * adding is not necessarily the right answer for
	     * concave polygons and subtractive triangles.
	     */
	    angle += M_PI ;
	} else if ( vec_dot < 1.0 ) {
	    if ( vec0[X] * vec1[Y] - vec1[X] * vec0[Y] >= 0.0 ) {
		/* add angle due to dot product of vectors */
		angle += acos( vec_dot ) ;
	    } else {
		/* subtract angle due to dot product of vectors */
		angle -= acos( vec_dot ) ;
	    }
	} /* if vec_dot >= 1.0, angle does not change */

	/* 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 */
    inside_flag = fmod( fabs(angle) + M_PI, 4.0*M_PI ) > 2.0*M_PI ;

    return( inside_flag ) ;
}

/* ======= Weiler algorithm ============================================ */

/* Track quadrant movements using Weiler's algorithm (elsewhere in Graphics
 * Gems IV).  Algorithm has been optimized for testing purposes, but the
 * crossings algorithm is still faster.	 Included to provide timings.
 *
 * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
 * _point_, returns 1 if inside, 0 if outside.	WINDING can be defined for
 * this test.
 */

#define QUADRANT( vtx, x, y )	\
	(((vtx)[X]>(x)) ? ( ((vtx)[Y]>(y)) ? 0:3 ) : ( ((vtx)[Y]>(y)) ? 1:2 ))

#define X_INTERCEPT( v0, v1, y )	\
	( (((v1)[X]-(v0)[X])/((v1)[Y]-(v0)[Y])) * ((y)-(v0)[Y]) + (v0)[X] )

int WeilerTest( pgon, numverts, point )
double	pgon[][2] ;
int	numverts ;
double	point[2] ;
{
register int	angle, qd, next_qd, delta, j ;
register double ty, tx, *vtx0, *vtx1 ;
int	inside_flag ;

    tx = point[X] ;
    ty = point[Y] ;

    vtx0 = pgon[numverts-1] ;
    qd = QUADRANT( vtx0, tx, ty ) ;
    angle = 0 ;

    vtx1 = pgon[0] ;

    for ( j = numverts+1 ; --j ; ) {
	/* calculate quadrant and delta from last quadrant */
	next_qd = QUADRANT( vtx1, tx, ty ) ;
	delta = next_qd - qd ;

	/* adjust delta and add it to total angle sum */
	switch ( delta ) {
	    case 0:	/* do nothing */
		break ;
	    case -1:
	    case 3:
		angle-- ;
		qd = next_qd ;
		break ;

	    case 1:
	    case -3:
		angle++ ;
		qd = next_qd ;
		break ;

	    case 2:
	    case -2:
		if (X_INTERCEPT( vtx0, vtx1, ty ) > tx ) {
		    angle -= delta ;
		} else {
		    angle += delta ;
		}
		qd = next_qd ;
		break ;
	}

	/* increment for next step */
	vtx0 = vtx1 ;
	vtx1 += 2 ;
    }

#ifdef	WINDING
    /* simple windings test:  if angle != 0, then point is inside */
    inside_flag = ( angle != 0 ) ;
#else
    /* Jordan test:  if angle is +-4, 12, 20, etc then point is inside */
    inside_flag = ( (angle/4) & 0x1 ) ;
#endif

    return (inside_flag) ;
}

#undef	QUADRANT
#undef	Y_INTERCEPT

/* ======= Triangle half-plane algorithm ================================== */

/* Split the polygon into a fan of triangles and for each triangle test if
 * the point is inside of the three half-planes formed by the triangle's edges.
 *
 * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
 * which returns a pointer to a plane set array.
 * Call testing procedure with a pointer to this array, _numverts_, and
 * test point _point_, returns 1 if inside, 0 if outside.
 * Call cleanup with pointer to plane set array to free space.
 *
 * SORT and CONVEX can be defined for this test.
 */

/* split polygons along set of x axes - call preprocess once */
pPlaneSet	PlaneSetup( pgon, numverts )
double	pgon[][2] ;
int	numverts ;
{
int	i, p1, p2 ;
double	tx, ty, vx0, vy0 ;
pPlaneSet	pps, pps_return ;
#ifdef	SORT
double	len[3], len_temp ;
int	j ;
PlaneSet	ps_temp ;
#ifdef	CONVEX
pPlaneSet	pps_new ;
pSizePlanePair p_size_pair ;
#endif
#endif

    pps = pps_return =
	    (pPlaneSet)malloc( 3 * (numverts-2) * sizeof( PlaneSet )) ;
    MALLOC_CHECK( pps ) ;
#ifdef	CONVEX
#ifdef	SORT
    p_size_pair =
	(pSizePlanePair)malloc( (numverts-2) * sizeof( SizePlanePair ) ) ;
    MALLOC_CHECK( p_size_pair ) ;
#endif
#endif

    vx0 = pgon[0][X] ;
    vy0 = pgon[0][Y] ;

    for ( p1 = 1, p2 = 2 ; p2 < numverts ; p1++, p2++ ) {
	pps->vx = vy0 - pgon[p1][Y] ;
	pps->vy = pgon[p1][X] - vx0 ;
	pps->c = pps->vx * vx0 + pps->vy * vy0 ;
#ifdef	SORT
	len[0] = pps->vx * pps->vx + pps->vy * pps->vy ;
#ifdef	CONVEX
#ifdef	HYBRID
	pps->ext_flag = ( p1 == 1 ) ;
#endif
	/* Sort triangles by areas, so compute (twice) the area here */
	p_size_pair[p1-1].pps = pps ;
	p_size_pair[p1-1].size =
			( pgon[0][X] * pgon[p1][Y] ) +
			( pgon[p1][X] * pgon[p2][Y] ) +
			( pgon[p2][X] * pgon[0][Y] ) -
			( pgon[p1][X] * pgon[0][Y] ) -
			( pgon[p2][X] * pgon[p1][Y] ) -
			( pgon[0][X] * pgon[p2][Y] ) ;
#endif
#endif
	pps++ ;
	pps->vx = pgon[p1][Y] - pgon[p2][Y] ;
	pps->vy = pgon[p2][X] - pgon[p1][X] ;
	pps->c = pps->vx * pgon[p1][X] + pps->vy * pgon[p1][Y] ;
#ifdef	SORT
	len[1] = pps->vx * pps->vx + pps->vy * pps->vy ;
#endif
#ifdef	CONVEX
#ifdef	HYBRID
	pps->ext_flag = TRUE ;
#endif
#endif
	pps++ ;
	pps->vx = pgon[p2][Y] - vy0 ;
	pps->vy = vx0 - pgon[p2][X] ;
	pps->c = pps->vx * pgon[p2][X] + pps->vy * pgon[p2][Y] ;
#ifdef	SORT
	len[2] = pps->vx * pps->vx + pps->vy * pps->vy ;
#endif
#ifdef	CONVEX
#ifdef	HYBRID
	pps->ext_flag = ( p2 == numverts-1 ) ;
#endif
#endif

	/* find an average point which must be inside of the triangle */
	tx = ( vx0 + pgon[p1][X] + pgon[p2][X] ) / 3.0 ;
	ty = ( vy0 + pgon[p1][Y] + pgon[p2][Y] ) / 3.0 ;

	/* check sense and reverse if test point is not thought to be inside
	 * first triangle
	 */
	if ( pps->vx * tx + pps->vy * ty >= pps->c ) {
	    /* back up to start of plane set */
	    pps -= 2 ;
	    /* point is thought to be outside, so reverse sense of edge
	     * normals so that it is correctly considered inside.
	     */
	    for ( i = 0 ; i < 3 ; i++ ) {
		pps->vx = -pps->vx ;
		pps->vy = -pps->vy ;
		pps->c	= -pps->c ;
		pps++ ;
	    }
	} else {
	    pps++ ;
	}

#ifdef	SORT
	/* sort the planes based on the edge lengths */
	pps -= 3 ;
	for ( i = 0 ; i < 2 ; i++ ) {
	    for ( j = i+1 ; j < 3 ; j++ ) {
		if ( len[i] < len[j] ) {
		    ps_temp = pps[i] ;
		    pps[i] = pps[j] ;
		    pps[j] = ps_temp ;
		    len_temp = len[i] ;
		    len[i] = len[j] ;
		    len[j] = len_temp ;
		}
	    }
	}
	pps += 3 ;
#endif
    }

#ifdef	CONVEX
#ifdef	SORT
    /* sort the triangles based on their areas */
    qsort( p_size_pair, numverts-2,
	    sizeof( SizePlanePair ), CompareSizePlanePairs ) ;

    /* make the plane sets match the sorted order */
    for ( i = 0, pps = pps_return
	; i < numverts-2
	; i++ ) {

	pps_new = p_size_pair[i].pps ;
	for ( j = 0 ; j < 3 ; j++, pps++, pps_new++ ) {
	    ps_temp = *pps ;
	    *pps = *pps_new ;
	    *pps_new = ps_temp ;
	}
    }
    free( p_size_pair ) ;
#endif
#endif

    return( pps_return ) ;
}

#ifdef	CONVEX
#ifdef	SORT
int CompareSizePlanePairs( p_sp0, p_sp1 )
pSizePlanePair	p_sp0, p_sp1 ;
{
    if ( p_sp0->size == p_sp1->size ) {
	return( 0 ) ;
    } else {
	return( p_sp0->size > p_sp1->size ? -1 : 1 ) ;
    }
}
#endif
#endif


/* check point for inside of three "planes" formed by triangle edges */
int PlaneTest( p_plane_set, numverts, point )
pPlaneSet	p_plane_set ;
int	numverts ;
double	point[2] ;
{
register pPlaneSet	ps ;
register int	p2 ;
#ifndef CONVEX
register int	inside_flag ;
#endif
register double tx, ty ;

    tx = point[X] ;
    ty = point[Y] ;

#ifndef CONVEX
    inside_flag = 0 ;
#endif

    for ( ps = p_plane_set, p2 = numverts-1 ; --p2 ; ) {

	if ( ps->vx * tx + ps->vy * ty < ps->c ) {
	    ps++ ;
	    if ( ps->vx * tx + ps->vy * ty < ps->c ) {
		ps++ ;
		/* note: we make the third edge have a slightly different
		 * equality condition, since this third edge is in fact
		 * the next triangle's first edge.  Not fool-proof, but
		 * it doesn't hurt (better would be to keep track of the
		 * triangle's area sign so we would know which kind of
		 * triangle this is).  Note that edge sorting nullifies
		 * this special inequality, too.
		 */
		if ( ps->vx * tx + ps->vy * ty <= ps->c ) {
		    /* point is inside polygon */
#ifdef CONVEX
		    return( 1 ) ;
#else
		    inside_flag = !inside_flag ;
#endif
		}
#ifdef	CONVEX
#ifdef	HYBRID
		/* check if outside exterior edge */
		else if ( ps->ext_flag ) return( 0 ) ;
#endif
#endif
		ps++ ;
	    } else {
#ifdef	CONVEX
#ifdef	HYBRID
		/* check if outside exterior edge */
		if ( ps->ext_flag ) return( 0 ) ;
#endif
#endif
		/* get past last two plane tests */
		ps += 2 ;