ptinpoly.c

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

	    }
	} else {
#ifdef	CONVEX
#ifdef	HYBRID
	    /* check if outside exterior edge */
	    if ( ps->ext_flag ) return( 0 ) ;
#endif
#endif
	    /* get past all three plane tests */
	    ps += 3 ;
	}
    }

#ifdef CONVEX
    /* for convex, if we make it to here, all triangles were missed */
    return( 0 ) ;
#else
    return( inside_flag ) ;
#endif
}

void PlaneCleanup( p_plane_set )
pPlaneSet	p_plane_set ;
{
    free( p_plane_set ) ;
}

/* ======= Barycentric algorithm ========================================== */

/* Split the polygon into a fan of triangles and for each triangle test if
 * the point has barycentric coordinates which are inside the triangle.
 * Similar to Badouel's code in Graphics Gems, with a little more efficient
 * coding.
 *
 * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
 * _point_, returns 1 if inside, 0 if outside.
 */
int BarycentricTest( pgon, numverts, point )
double	pgon[][2] ;
int	numverts ;
double	point[2] ;
{
register double *pg1, *pg2, *pgend ;
register double tx, ty, u0, u1, u2, v0, v1, vx0, vy0, alpha, beta, denom ;
int	inside_flag ;

    tx = point[X] ;
    ty = point[Y] ;
    vx0 = pgon[0][X] ;
    vy0 = pgon[0][Y] ;
    u0 = tx - vx0 ;
    v0 = ty - vy0 ;

    inside_flag = 0 ;
    pgend = pgon[numverts-1] ;
    for ( pg1 = pgon[1], pg2 = pgon[2] ; pg1 != pgend ; pg1+=2, pg2+=2 ) {

	u1 = pg1[X] - vx0 ;
	if ( u1 == 0.0 ) {

	    /* 0 and 1 vertices have same X value */

	    /* zero area test - can be removed for convex testing */
	    u2 = pg2[X] - vx0 ;
	    if ( ( u2 == 0.0 ) ||

		/* compute beta and check bounds */
		/* we use "<= 0.0" so that points on the shared interior
		 * edge will (generally) be inside only one polygon.
		 */
		 ( ( beta = u0 / u2 ) <= 0.0 ) ||
		 ( beta > 1.0 ) ||

		 /* zero area test - remove for convex testing */
		 ( ( v1 = pg1[Y] - vy0 ) == 0.0 ) ||

		 /* compute alpha and check bounds */
		 ( ( alpha = ( v0 - beta *
		    ( pg2[Y] - vy0 ) ) / v1 ) < 0.0 ) ) {

		/* whew! missed! */
		goto NextTri ;
	    }

	} else {
	    /* 0 and 1 vertices have different X value */

	    /* compute denom and check for zero area triangle - check
	     * is not needed for convex polygon testing
	     */
	    u2 = pg2[X] - vx0 ;
	    v1 = pg1[Y] - vy0 ;
	    denom = ( pg2[Y] - vy0 ) * u1 - u2 * v1 ;
	    if ( ( denom == 0.0 ) ||

		/* compute beta and check bounds */
		/* we use "<= 0.0" so that points on the shared interior
		 * edge will (generally) be inside only one polygon.
		 */
		 ( ( beta = ( v0 * u1 - u0 * v1 ) / denom ) <= 0.0 ) ||
		 ( beta > 1.0 ) ||

		 /* compute alpha & check bounds */
		 ( ( alpha = ( u0 - beta * u2 ) / u1 ) < 0.0 ) ) {

		/* whew! missed! */
		goto NextTri ;
	    }
	}

	/* check gamma */
	if ( alpha + beta <= 1.0 ) {
	    /* survived */
	    inside_flag = !inside_flag ;
	}

	NextTri: ;
    }
    return( inside_flag ) ;
}

/* ======= Barycentric precompute (Spackman) algorithm ===================== */

/* Split the polygon into a fan of triangles and for each triangle test if
 * the point has barycentric coordinates which are inside the triangle.
 * Use Spackman's normalization method to precompute various parameters.
 *
 * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
 * which returns a pointer to the array of the parameters records and the
 * number of parameter records created.
 * Call testing procedure with the first vertex in the polygon _pgon[0]_,
 * a pointer to this array, the number of parameter records, and test point
 * _point_, returns 1 if inside, 0 if outside.
 * Call cleanup with pointer to parameter record array to free space.
 *
 * SORT can be defined for this test.
 * (CONVEX could be added: see PlaneSetup and PlaneTest for method)
 */
pSpackmanSet	SpackmanSetup( pgon, numverts, p_numrec )
double	pgon[][2] ;
int	numverts ;
int	*p_numrec ;
{
int	p1, p2, degen ;
double	denom, u1, v1, *pv[3] ;
pSpackmanSet	pss, pss_return ;
#ifdef	SORT
double	u[2], v[2], len[2], *pv_temp ;
#endif

    pss = pss_return =
	    (pSpackmanSet)malloc( (numverts-2) * sizeof( SpackmanSet )) ;
    MALLOC_CHECK( pss ) ;

    degen = 0 ;

    for ( p1 = 1, p2 = 2 ; p2 < numverts ; p1++, p2++ ) {

	pv[0] = pgon[0] ;
	pv[1] = pgon[p1] ;
	pv[2] = pgon[p2] ;

#ifdef	SORT
	/* Note that sorting can cause a mismatch of alpha/beta inequality
	 * tests.  In other words, test points on an interior line between
	 * test triangles will often then be wrong.
	 */
	u[0] = pv[1][X] - pv[0][X] ;
	u[1] = pv[2][X] - pv[0][X] ;
	v[0] = pv[1][Y] - pv[0][Y] ;
	v[1] = pv[2][Y] - pv[0][Y] ;
	len[0] = u[0] * u[0] + v[0] * v[0] ;
	len[1] = u[1] * u[1] + v[1] * v[1] ;

	/* compare two edges touching anchor point and put longest first */
	/* we don't sort all three edges because the anchor point and
	 * values computed from it gets used for all triangles in the fan.
	 */
	if ( len[0] < len[1] ) {
	    pv_temp = pv[1] ;
	    pv[1] = pv[2] ;
	    pv[2] = pv_temp ;
	}
#endif

	u1 = pv[1][X] - pv[0][X] ;
	pss->u2 = pv[2][X] - pv[0][X] ;
	v1 = pv[1][Y] - pv[0][Y] ;
	pss->v2 = pv[2][Y] - pv[0][Y] ;
	pss->u1_nonzero = !( u1 == 0.0 ) ;
	if ( pss->u1_nonzero ) {
	    /* not zero, so compute inverse */
	    pss->inv_u1 = 1.0 / u1 ;
	    denom = pss->v2 * u1 - pss->u2 * v1 ;
	    if ( denom == 0.0 ) {
		/* degenerate triangle, ignore it */
		degen++ ;
		goto Skip ;
	    } else {
		pss->u1p = u1 / denom ;
		pss->v1p = v1 / denom ;
	    }
	} else {
	    if ( pss->u2 == 0.0 ) {
		/* degenerate triangle, ignore it */
		degen++ ;
		goto Skip ;
	    } else {
		/* not zero, so compute inverse */
		pss->inv_u2 = 1.0 / pss->u2 ;
		if ( v1 == 0.0 ) {
		    /* degenerate triangle, ignore it */
		    degen++ ;
		    goto Skip ;
		} else {
		    pss->inv_v1 = 1.0 / v1 ;
		}
	    }
	}

	pss++ ;
	Skip: ;
    }

    /* number of Spackman records */
    *p_numrec = numverts - degen - 2 ;
    if ( degen ) {
	pss = pss_return =
		(pSpackmanSet)realloc( pss_return,
			(numverts-2-degen) * sizeof( SpackmanSet )) ;
    }

    return( pss_return ) ;
}

/* barycentric, a la Gems I and Spackman's normalization precompute */
int SpackmanTest( anchor, p_spackman_set, numrec, point )
double	anchor[2] ;
pSpackmanSet	p_spackman_set ;
int	numrec ;
double	point[2] ;
{
register pSpackmanSet	pss ;
register int	inside_flag ;
register int	nr ;
register double tx, ty, vx0, vy0, u0, v0, alpha, beta ;

    tx = point[X] ;
    ty = point[Y] ;
    /* note that we really need only the first vertex of the polygon,
     * so do not really need to keep the whole polygon around.
     */
    vx0 = anchor[X] ;
    vy0 = anchor[Y] ;
    u0 = tx - vx0 ;
    v0 = ty - vy0 ;

    inside_flag = 0 ;

    for ( pss = p_spackman_set, nr = numrec+1 ; --nr ; pss++ ) {

	if ( pss->u1_nonzero ) {
	    /* 0 and 2 vertices have different X value */

	    /* compute beta and check bounds */
	    /* we use "<= 0.0" so that points on the shared interior edge
	     * will (generally) be inside only one polygon.
	     */
	    beta = ( v0 * pss->u1p - u0 * pss->v1p ) ;
	    if ( ( beta <= 0.0 ) || ( beta > 1.0 ) ||

		 /* compute alpha & check bounds */
		 ( ( alpha = ( u0 - beta * pss->u2 ) * pss->inv_u1 )
			< 0.0 ) ) {

		/* whew! missed! */
		goto NextTri ;
	    }
	} else {
	    /* 0 and 2 vertices have same X value */

	    /* compute beta and check bounds */
	    /* we use "<= 0.0" so that points on the shared interior edge
	     * will (generally) be inside only one polygon.
	     */
	    beta = u0 * pss->inv_u2 ;
	    if ( ( beta <= 0.0 ) || ( beta >= 1.0 ) ||

		 /* compute alpha and check bounds */
		 ( ( alpha = ( v0 - beta * pss->v2 ) * pss->inv_v1 )
			< 0.0 ) ) {

		/* whew! missed! */
		goto NextTri ;
	    }
	}

	/* check gamma */
	if ( alpha + beta <= 1.0 ) {
	    /* survived */
	    inside_flag = !inside_flag ;
	}

	NextTri: ;
    }

    return( inside_flag ) ;
}

void SpackmanCleanup( p_spackman_set )
pSpackmanSet	p_spackman_set ;
{
    free( p_spackman_set ) ;
}

/* ======= Trapezoid (bin) algorithm ======================================= */

/* Split polygons along set of y bins and sorts the edge fragments.  Testing
 * is done against these fragments.
 *
 * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices, the
 * number of bins desired _bins_, and a pointer to a trapezoid structure
 * _p_trap_set_.
 * Call testing procedure with 2D polygon _pgon_ with _numverts_ number of
 * vertices, _p_trap_set_ pointer to trapezoid structure, and test point
 * _point_, returns 1 if inside, 0 if outside.
 * Call cleanup with pointer to trapezoid structure to free space.
 */
void TrapezoidSetup( pgon, numverts, bins, p_trap_set )
double	pgon[][2] ;
int	numverts ;
int	bins ;
pTrapezoidSet	p_trap_set ;
{
double	*vtx0, *vtx1, *vtxa, *vtxb, slope ;
int	i, j, bin_tot[TOT_BINS], ba, bb, id, full_cross, count ;
double	fba, fbb, vx0, vx1, dy, vy0 ;

    p_trap_set->bins = bins ;
    p_trap_set->trapz = (pTrapezoid)malloc( p_trap_set->bins *
	    sizeof(Trapezoid));
    MALLOC_CHECK( p_trap_set->trapz ) ;

    p_trap_set->minx =
    p_trap_set->maxx = pgon[0][X] ;
    p_trap_set->miny =
    p_trap_set->maxy = pgon[0][Y] ;

    for ( i = 1 ; i < numverts ; i++ ) {
	if ( p_trap_set->minx > (vx0 = pgon[i][X]) ) {
	    p_trap_set->minx = vx0 ;
	} else if ( p_trap_set->maxx < vx0 ) {
	    p_trap_set->maxx = vx0 ;
	}

	if ( p_trap_set->miny > (vy0 = pgon[i][Y]) ) {
	    p_trap_set->miny = vy0 ;
	} else if ( p_trap_set->maxy < vy0 ) {
	    p_trap_set->maxy = vy0 ;
	}
    }

    /* add a little to the bounds to ensure everything falls inside area */
    p_trap_set->miny -= EPSILON * (p_trap_set->maxy-p_trap_set->miny) ;
    p_trap_set->maxy += EPSILON * (p_trap_set->maxy-p_trap_set->miny) ;

    p_trap_set->ydelta =
	    (p_trap_set->maxy-p_trap_set->miny) / (double)p_trap_set->bins ;
    p_trap_set->inv_ydelta = 1.0 / p_trap_set->ydelta ;

    /* find how many locations to allocate for each bin */
    for ( i = 0 ; i < p_trap_set->bins ; i++ ) {
	bin_tot[i] = 0 ;
    }

    vtx0 = pgon[numverts-1] ;
    for ( i = 0 ; i < numverts ; i++ ) {
	vtx1 = pgon[i] ;

	/* skip if Y's identical (edge has no effect) */
	if ( vtx0[Y] != vtx1[Y] ) {

	    if ( vtx0[Y] < vtx1[Y] ) {
		vtxa = vtx0 ;
		vtxb = vtx1 ;
	    } else {
		vtxa = vtx1 ;
		vtxb = vtx0 ;
	    }
	    ba = (int)(( vtxa[Y]-p_trap_set->miny ) * p_trap_set->inv_ydelta) ;
	    fbb = ( vtxb[Y] - p_trap_set->miny ) * p_trap_set->inv_ydelta ;
	    bb = (int)fbb ;
	    /* if high vertex ends on a boundary, don't go into next boundary */
	    if ( fbb - (double)bb == 0.0 ) {
		bb-- ;
	    }

	    /* mark the bins with this edge */
	    for ( j = ba ; j <= bb ; j++ ) {
		bin_tot[j]++ ;
	    }
	}

	vtx0 = vtx1 ;
    }

    /* allocate the bin contents and fill in some basics */
    for ( i = 0 ; i < p_trap_set->bins ; i++ ) {
	p_trap_set->trapz[i].edge_set =
		(pEdge*)malloc( bin_tot[i] * sizeof(pEdge) ) ;
	MALLOC_CHECK( p_trap_set->trapz[i].edge_set ) ;
	for ( j = 0 ; j < bin_tot[i] ; j++ ) {
	    p_trap_set->trapz[i].edge_set[j] =
		(pEdge)malloc( sizeof(Edge) ) ;
	    MALLOC_CHECK( p_trap_set->trapz[i].edge_set[j] ) ;
	}

	/* start these off at some awful values; refined below */
	p_trap_set->trapz[i].minx = p_trap_set->maxx ;
	p_trap_set->trapz[i].maxx = p_trap_set->minx ;
	p_trap_set->trapz[i].count = 0 ;
    }

    /* now go through list yet again, putting edges in bins */
    vtx0 = pgon[numverts-1] ;
    id = numverts-1 ;
    for ( i = 0 ; i < numverts ; i++ ) {
	vtx1 = pgon[i] ;

	/* we can skip edge if Y's are equal */
	if ( vtx0[Y] != vtx1[Y] ) {
	    if ( vtx0[Y] < vtx1[Y] ) {
		vtxa = vtx0 ;
		vtxb = vtx1 ;
	    } else {
		vtxa = vtx1 ;
		vtxb = vtx0 ;
	    }
	    fba = ( vtxa[Y] - p_trap_set->miny ) * p_trap_set->inv_ydelta ;
	    ba = (int)fba ;
	    fbb = ( vtxb[Y] - p_trap_set->miny ) * p_trap_set->inv_ydelta ;
	    bb = (int)fbb ;
	    /* if high vertex ends on a boundary, don't go into it */
	    if ( fbb == (double)bb ) {
		bb-- ;
	    }

	    vx0 = vtxa[X] ;
	    dy = vtxa[Y] - vtxb[Y] ;
	    slope = p_trap_set->ydelta * ( vtxa[X] - vtxb[X] ) / dy ;

	    /* set vx1 in case loop is not entered */
	    vx1 = vx0 ;
	    full_cross = 0 ;

	    for ( j = ba ; j < bb ; j++, vx0 = vx1 ) {
		/* could increment vx1, but for greater accuracy recompute it */
		vx1 = vtxa[X] + ( (double)(j+1) - fba ) * slope ;

		count = p_trap_set->trapz[j].count++ ;
		p_trap_set->trapz[j].edge_set[count]->id = id ;
		p_trap_set->trapz[j].edge_set[count]->full_cross = full_cross;
		TrapBound( j, count, vx0, vx1, p_trap_set ) ;
		full_cross = 1 ;
	    }

	    /* at last bin - fill as above, but with vx1 = vtxb[X] */
	    vx0 = vx1 ;
	    vx1 = vtxb[X] ;
	    count = p_trap_set->trapz[bb].count++ ;
	    p_trap_set->trapz[bb].edge_set[count]->id = id ;
	    /* the last bin is never a full crossing */
	    p_trap_set->trapz[bb].edge_set[count]->full_cross = 0 ;
	    TrapBound( bb, count, vx0, vx1, p_trap_set ) ;
	}

	vtx0 = vtx1 ;
	id = i ;
    }

    /* finally, sort the bins' contents by minx */
    for ( i = 0 ; i < p_trap_set->bins ; i++ ) {
	qsort( p_trap_set->trapz[i].edge_set, p_trap_set->trapz[i].count,
		sizeof(pEdge), CompareEdges ) ;
    }
}

void TrapBound( j, count, vx0, vx1, p_trap_set )
int	j, count ;
double	vx0, vx1 ;
pTrapezoidSet	p_trap_set ;
{
double	xt ;

    if ( vx0 > vx1 ) {
	xt = vx0 ;
	vx0 = vx1 ;
	vx1 = xt ;
    }

    if ( p_trap_set->trapz[j].minx > vx0 ) {
	p_trap_set->trapz[j].minx = vx0 ;
    }
    if ( p_trap_set->trapz[j].maxx < vx1 ) {
	p_trap_set->trapz[j].maxx = vx1 ;
    }
    p_trap_set->trapz[j].edge_set[count]->minx = vx0 ;
    p_trap_set->trapz[j].edge_set[count]->maxx = vx1 ;
}

/* used by qsort to sort */
int CompareEdges( u, v )
pEdge	*u, *v ;
{
    if ( (*u)->minx == (*v)->minx ) {
	return( 0 ) ;
    } else {
	return( (*u)->minx < (*v)->minx ? -1 : 1 ) ;
    }
}

int TrapezoidTest( pgon, numverts, p_trap_set, point )
double	pgon[][2] ;
int	numverts ;
pTrapezoidSet	p_trap_set ;
double	point[2] ;
{
int	j, b, count, id ;
double	tx, ty, *vtx0, *vtx1 ;
pEdge	*pp_bin ;
pTrapezoid	p_trap ;
int	inside_flag ;

    inside_flag = 0 ;