ptinpoly.c

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

		for ( j = count+1 ; --j ; p_gr++ ) {
		    /* test if y is between edges */
		    if ( ty >= p_gr->miny && ty < p_gr->maxy ) {
			if ( tx > p_gr->maxx ) {
			    inside_flag = !inside_flag ;
			} else if ( tx > p_gr->minx ) {
			    /* full computation */
			    if ( ( p_gr->xa -
				( p_gr->ya - ty ) * p_gr->slope ) < tx ) {
				inside_flag = !inside_flag ;
			    }
			}
		    }
		}
		break ;

	    case GC_AIM_B:
		/* bottom edge is clear, shoot Y+ ray */
		/* note - this next statement requires that GC_BL_IN is 1 */
		inside_flag = gc_flags & GC_BL_IN ;
		for ( j = count+1 ; --j ; p_gr++ ) {
		    /* test if x is between edges */
		    if ( tx >= p_gr->minx && tx < p_gr->maxx ) {
			if ( ty > p_gr->maxy ) {
			    inside_flag = !inside_flag ;
			} else if ( ty > p_gr->miny ) {
			    /* full computation */
			    if ( ( p_gr->ya - ( p_gr->xa - tx ) *
				    p_gr->inv_slope ) < ty ) {
				inside_flag = !inside_flag ;
			    }
			}
		    }
		}
		break ;

	    case GC_AIM_R:
		/* right edge is clear, shoot X+ ray */
		inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0 ;

		/* TBD: Note, we could have sorted the edges to be tested
		 * by miny or somesuch, and so be able to cut testing
		 * short when the list's miny > point.y .
		 */
		for ( j = count+1 ; --j ; p_gr++ ) {
		    /* test if y is between edges */
		    if ( ty >= p_gr->miny && ty < p_gr->maxy ) {
			if ( tx <= p_gr->minx ) {
			    inside_flag = !inside_flag ;
			} else if ( tx <= p_gr->maxx ) {
			    /* full computation */
			    if ( ( p_gr->xa -
				( p_gr->ya - ty ) * p_gr->slope ) >= tx ) {
				inside_flag = !inside_flag ;
			    }
			}
		    }
		}
		break ;

	    case GC_AIM_T:
		/* top edge is clear, shoot Y+ ray */
		inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0 ;
		for ( j = count+1 ; --j ; p_gr++ ) {
		    /* test if x is between edges */
		    if ( tx >= p_gr->minx && tx < p_gr->maxx ) {
			if ( ty <= p_gr->miny ) {
			    inside_flag = !inside_flag ;
			} else if ( ty <= p_gr->maxy ) {
			    /* full computation */
			    if ( ( p_gr->ya - ( p_gr->xa - tx ) *
				    p_gr->inv_slope ) >= ty ) {
				inside_flag = !inside_flag ;
			    }
			}
		    }
		}
		break ;

	    case GC_AIM_C:
		/* no edge is clear, bite the bullet and test
		 * against the bottom left corner.
		 * We use Franklin Antonio's algorithm (Graphics Gems III).
		 */
		/* TBD: Faster yet might be to test against the closest
		 * corner to the cell location, but our hope is that we
		 * rarely need to do this testing at all.
		 */
		inside_flag = ((gc_flags & GC_BL_IN) == GC_BL_IN) ;
		init_flag = TRUE ;

		/* get lower left corner coordinate */
		cornerx = p_gs->glx[(int)xcell] ;
		cornery = p_gs->gly[(int)ycell] ;
		for ( j = count+1 ; --j ; p_gr++ ) {

		    /* quick out test: if test point is
		     * less than minx & miny, edge cannot overlap.
		     */
		    if ( tx >= p_gr->minx && ty >= p_gr->miny ) {

			/* quick test failed, now check if test point and
			 * corner are on different sides of edge.
			 */
			if ( init_flag ) {
			    /* Compute these at most once for test */
			    /* P3 - P4 */
			    bx = tx - cornerx ;
			    by = ty - cornery ;
			    init_flag = FALSE ;
			}
			denom = p_gr->ay * bx - p_gr->ax * by ;
			if ( denom != 0.0 ) {
			    /* lines are not collinear, so continue */
			    /* P1 - P3 */
			    cx = p_gr->xa - tx ;
			    cy = p_gr->ya - ty ;
			    alpha = by * cx - bx * cy ;
			    if ( denom > 0.0 ) {
				if ( alpha < 0.0 || alpha >= denom ) {
				    /* test edge not hit */
				    goto NextEdge ;
				}
				beta = p_gr->ax * cy - p_gr->ay * cx ;
				if ( beta < 0.0 || beta >= denom ) {
				    /* polygon edge not hit */
				    goto NextEdge ;
				}
			    } else {
				if ( alpha > 0.0 || alpha <= denom ) {
				    /* test edge not hit */
				    goto NextEdge ;
				}
				beta = p_gr->ax * cy - p_gr->ay * cx ;
				if ( beta > 0.0 || beta <= denom ) {
				    /* polygon edge not hit */
				    goto NextEdge ;
				}
			    }
			    inside_flag = !inside_flag ;
			}

		    }
		    NextEdge: ;
		}
		break ;
	    }
	} else {
	    /* simple cell, so if lower left corner is in,
	     * then cell is inside.
	     */
	    inside_flag = p_gc->gc_flags & GC_BL_IN ;
	}
    }

    return( inside_flag ) ;
}

void GridCleanup( p_gs )
pGridSet	p_gs ;
{
int	i ;
pGridCell	p_gc ;

    for ( i = 0, p_gc = p_gs->gc
	; i < p_gs->tot_cells
	; i++, p_gc++ ) {

	if ( p_gc->gr ) {
	    free( p_gc->gr ) ;
	}
    }
    free( p_gs->glx ) ;
    free( p_gs->gly ) ;
    free( p_gs->gc ) ;
}

/* ======= Exterior (convex only) algorithm =============================== */

/* Test the edges of the convex polygon against the point.  If the point is
 * outside any edge, the point is outside the polygon.
 *
 * 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.
 *
 * RANDOM can be defined for this test.
 * CONVEX must be defined for this test; it is not usable for general polygons.
 */

#ifdef	CONVEX
/* make exterior plane set */
pPlaneSet ExteriorSetup( pgon, numverts )
double	pgon[][2] ;
int	numverts ;
{
int	p1, p2, flip_edge ;
pPlaneSet	pps, pps_return ;
#ifdef	RANDOM
int	i, ind ;
PlaneSet	ps_temp ;
#endif

    pps = pps_return =
	    (pPlaneSet)malloc( numverts * sizeof( PlaneSet )) ;
    MALLOC_CHECK( pps ) ;

    /* take cross product of vertex to find handedness */
    flip_edge = (pgon[0][X] - pgon[1][X]) * (pgon[1][Y] - pgon[2][Y] ) >
		(pgon[0][Y] - pgon[1][Y]) * (pgon[1][X] - pgon[2][X] ) ;

    /* Generate half-plane boundary equations now for faster testing later.
     * vx & vy are the edge's normal, c is the offset from the origin.
     */
    for ( p1 = numverts-1, p2 = 0 ; p2 < numverts ; p1 = p2, p2++, 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] ;

	/* check sense and reverse plane edge if need be */
	if ( flip_edge ) {
	    pps->vx = -pps->vx ;
	    pps->vy = -pps->vy ;
	    pps->c  = -pps->c ;
	}
    }

#ifdef	RANDOM
    /* Randomize the order of the edges to improve chance of early out */
    /* There are better orders, but the default order is the worst */
    for ( i = 0, pps = pps_return
	; i < numverts
	; i++ ) {

	ind = (int)(RAN01() * numverts ) ;
	if ( ( ind < 0 ) || ( ind >= numverts ) ) {
	    fprintf( stderr,
		"Yikes, the random number generator is returning values\n" ) ;
	    fprintf( stderr,
		"outside the range [0.0,1.0), so please fix the code!\n" ) ;
	    ind = 0 ;
	}

	/* swap edges */
	ps_temp = *pps ;
	*pps = pps_return[ind] ;
	pps_return[ind] = ps_temp ;
    }
#endif

    return( pps_return ) ;
}

/* Check point for outside of all planes */
/* note that we don't need "pgon", since it's been processed into
 * its corresponding PlaneSet.
 */
int ExteriorTest( p_ext_set, numverts, point )
pPlaneSet	p_ext_set ;
int	numverts ;
double	point[2] ;
{
register PlaneSet	*pps ;
register int	p0 ;
register double tx, ty ;
int	inside_flag ;

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

    for ( p0 = numverts+1, pps = p_ext_set ; --p0 ; pps++ ) {

	/* test if the point is outside this edge */
	if ( pps->vx * tx + pps->vy * ty > pps->c ) {
	    return( 0 ) ;
	}
    }
    /* if we make it to here, we were inside all edges */
    return( 1 ) ;
}

void ExteriorCleanup( p_ext_set )
pPlaneSet	p_ext_set ;
{
    free( p_ext_set ) ;
}
#endif

/* ======= Inclusion (convex only) algorithm ============================== */

/* Create an efficiency structure (see Preparata) for the convex polygon which
 * allows binary searching to find which edge to test the point against.  This
 * algorithm is O(log n).
 *
 * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
 * which returns a pointer to an inclusion anchor structure.
 * Call testing procedure with a pointer to this structure and test point
 * _point_, returns 1 if inside, 0 if outside.
 * Call cleanup with pointer to inclusion anchor structure to free space.
 *
 * CONVEX must be defined for this test; it is not usable for general polygons.
 */

#ifdef	CONVEX
/* make inclusion wedge set */
pInclusionAnchor InclusionSetup( pgon, numverts )
double	pgon[][2] ;
int	numverts ;
{
int	pc, p1, p2, flip_edge ;
double	ax,ay, qx,qy, wx,wy, len ;
pInclusionAnchor pia ;
pInclusionSet	pis ;

    /* double the first edge to avoid needing modulo during test search */
    pia = (pInclusionAnchor)malloc( sizeof( InclusionAnchor )) ;
    MALLOC_CHECK( pia ) ;
    pis = pia->pis =
	    (pInclusionSet)malloc( (numverts+1) * sizeof( InclusionSet )) ;
    MALLOC_CHECK( pis ) ;

    pia->hi_start = numverts - 1 ;

    /* get average point to make wedges from */
    qx = qy = 0.0 ;
    for ( p2 = 0 ; p2 < numverts ; p2++ ) {
	qx += pgon[p2][X] ;
	qy += pgon[p2][Y] ;
    }
    pia->qx = qx /= (double)numverts ;
    pia->qy = qy /= (double)numverts ;

    /* take cross product of vertex to find handedness */
    pia->flip_edge = flip_edge =
		(pgon[0][X] - pgon[1][X]) * (pgon[1][Y] - pgon[2][Y] ) >
		(pgon[0][Y] - pgon[1][Y]) * (pgon[1][X] - pgon[2][X] ) ;


    ax = pgon[0][X] - qx ;
    ay = pgon[0][Y] - qy ;
    len = sqrt( ax * ax + ay * ay ) ;
    if ( len == 0.0 ) {
	fprintf( stderr, "sorry, polygon for inclusion test is defective\n" ) ;
	exit(1) ;
    }
    pia->ax = ax /= len ;
    pia->ay = ay /= len ;

    /* loop through edges, and double last edge */
    for ( pc = p1 = 0, p2 = 1
	; pc <= numverts
	; pc++, p1 = p2, p2 = (++p2)%numverts, pis++ ) {

	/* wedge border */
	wx = pgon[p1][X] - qx ;
	wy = pgon[p1][Y] - qy ;
	len = sqrt( wx * wx + wy * wy ) ;
	wx /= len ;
	wy /= len ;

	/* cosine of angle from anchor border to wedge border */
	pis->dot = ax * wx + ay * wy ;
	/* sign from cross product */
	if ( ( ax * wy > ay * wx ) == flip_edge ) {
	    pis->dot = -2.0 - pis->dot ;
	}

	/* edge */
	pis->ex = pgon[p1][Y] - pgon[p2][Y] ;
	pis->ey = pgon[p2][X] - pgon[p1][X] ;
	pis->ec = pis->ex * pgon[p1][X] + pis->ey * pgon[p1][Y] ;

	/* check sense and reverse plane eqns if need be */
	if ( flip_edge ) {
	    pis->ex = -pis->ex ;
	    pis->ey = -pis->ey ;
	    pis->ec = -pis->ec ;
	}
    }
    /* set first angle a little > 1.0 and last < -3.0 just to be safe. */
    pia->pis[0].dot = -3.001 ;
    pia->pis[numverts].dot = 1.001 ;

    return( pia ) ;
}

/* Find wedge point is in by binary search, then test wedge */
int InclusionTest( pia, point )
pInclusionAnchor	pia ;
double	point[2] ;
{
register double tx, ty, len, dot ;
int	inside_flag, lo, hi, ind ;
pInclusionSet	pis ;

    tx = point[X] - pia->qx ;
    ty = point[Y] - pia->qy ;
    len = sqrt( tx * tx + ty * ty ) ;
    /* check if point is exactly at anchor point (which is inside polygon) */
    if ( len == 0.0 ) return( 1 ) ;
    tx /= len ;
    ty /= len ;

    /* get dot product for searching */
    dot = pia->ax * tx + pia->ay * ty ;
    if ( ( pia->ax * ty > pia->ay * tx ) == pia->flip_edge ) {
	dot = -2.0 - dot ;
    }

    /* binary search through angle list and find matching angle pair */
    lo = 0 ;
    hi = pia->hi_start ;
    while ( lo <= hi ) {
	ind = (lo+hi)/ 2 ;
	if ( dot < pia->pis[ind].dot ) {
	    hi = ind - 1 ;
	} else if ( dot > pia->pis[ind+1].dot ) {
	    lo = ind + 1 ;
	} else {
	    goto Foundit ;
	}
    }
    /* should never reach here, but just in case... */
    fprintf( stderr,
	    "Hmmm, something weird happened - bad dot product %lg\n", dot);

    Foundit:

    /* test if the point is outside the wedge's exterior edge */
    pis = &pia->pis[ind] ;
    inside_flag = ( pis->ex * point[X] + pis->ey * point[Y] <= pis->ec ) ;

    return( inside_flag ) ;
}

void InclusionCleanup( p_inc_anchor )
pInclusionAnchor p_inc_anchor ;
{
    free( p_inc_anchor->pis ) ;
    free( p_inc_anchor ) ;
}
#endif


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

/*
 * This version is usually somewhat faster than the original published in
 * Graphics Gems IV; by turning the division for testing the X axis crossing
 * into a tricky multiplication test this part of the test became faster,
 * which had the additional effect of making the test for "both to left or
 * both to right" a bit slower for triangles than simply computing the
 * intersection each time.  The main increase is in triangle testing speed,
 * which was about 15% faster; all other polygon complexities were pretty much
 * the same as before.  On machines where division is very expensive (not the
 * case on the HP 9000 series on which I tested) this test should be much
 * faster overall than the old code.  Your mileage may (in fact, will) vary,
 * depending on the machine and the test data, but in general I believe this
 * code is both shorter and faster.  This test was inspired by unpublished
 * Graphics Gems submitted by Joseph Samosky and Mark Haigh-Hutchinson.
 * Related work by Samosky is in:
 *
 * Samosky, Joseph, "SectionView: A system for interactively specifying and
 * visualizing sections through three-dimensional medical image data",
 * M.S. Thesis, Department of Electrical Engineering and Computer Science,
 * Massachusetts Institute of Technology, 1993.
 *
 */

/* Shoot a test ray along +X axis.  The strategy 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.  Note that CONVEX and WINDING code can be added as
 * for the CrossingsTest() code; it is left out here for clarity.
 *
 * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
 * _point_, returns 1 if inside, 0 if outside.
 */
int CrossingsMultiplyTest( pgon, numverts, point )
double	pgon[][2] ;
int	numverts ;
double	point[2] ;
{
register int	j, yflag0, yflag1, inside_flag ;
register double	ty, tx, *vtx0, *vtx1 ;

    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] ;

    inside_flag = 0 ;
    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.
	 * The old test also checked whether the endpoints are both to the
	 * right or to the left of the test point.  However, given the faster
	 * intersection point computation used below, this test was found to
	 * be a break-even proposition for most polygons and a loser for
	 * triangles (where 50% or more of the edges which survive this test
	 * will cross quadrants and so have to have the X intersection computed
	 * anyway).  I credit Joseph Samosky with inspiring me to try dropping
	 * the "both left or both right" part of my code.
	 */
	if ( yflag0 != yflag1 ) {
	    /* Check intersection of pgon segment with +X ray.
	     * Note if >= point's X; if so, the ray hits it.
	     * The division operation is avoided for the ">=" test by checking
	     * the sign of the first vertex wrto the test point; idea inspired
	     * by Joseph Samosky's and Mark Haigh-Hutchinson's different
	     * polygon inclusion tests.
	     */
	    if ( ((vtx1[Y]-ty) * (vtx0[X]-vtx1[X]) >=
		    (vtx1[X]-tx) * (vtx0[Y]-vtx1[Y])) == yflag1 ) {
		inside_flag = !inside_flag ;
	    }
	}

	/* Move to the next pair of vertices, retaining info as possible. */
	yflag0 = yflag1 ;
	vtx0 = vtx1 ;
	vtx1 += 2 ;
    }

    return( inside_flag ) ;
}