p_test.c

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

/* p_test.c - point in polygon inside/outside tester.  This is
 * simply testing and display code, the actual algorithms are in ptinpoly.c.
 *
 * Probably the most important thing to set for timings is TEST_RATIO (too
 * low and the timings are untrustworthy, too high and you wait forever).
 * Start low and see how consistent separate runs appear to be.
 *
 * To add a new algorithm to the test suite, add code at the spots marked
 * with '+++'.
 *
 * See Usage() for command line options (or just do "p_test -?").
 *
 * by Eric Haines, 3D/Eye Inc, [e-mail]
 */

/* Define TIMER to perform timings on code (need system timing function) */
/* #define TIMER */

/* Number of times to try a single point vs. a polygon, per vertex.
 * This should be greater than 1 / ( HZ * approx. single test time in seconds )
 * in order to get a meaningful timings difference.  200 is reasonable for a
 * IBM PC 25 MHz 386 with no FPU, 50000 is good for an HP 720 workstation.
 * Start low and see how consistent separate runs appear to be.
 */

#ifndef M_PI
#define M_PI	3.14159265358979323846
#endif

#ifdef	TIMER
#define MACHINE_TEST_RATIO	50000
#else
#define MACHINE_TEST_RATIO	    1
#endif

/* ===========	that's all the easy stuff than can be changed  ============ */

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

#ifdef	TIMER
#ifdef	IBM_PC
#include <bios.h>
#include <time.h>
#define HZ	CLK_TCK
#define mGetTime(t)	(t) = biostime(0,0L) ;

#else	/* UNIX */
#include <sys/types.h>
#include <sys/param.h>
#include <sys/times.h>
struct	tms	Timebuf ;

#define mGetTime(t)	(t) = times(&Timebuf) ;
#if (!(defined(sparc)||defined(__sparc__)))
#define mGetTime(t)     (t) = times(&Timebuf) ;
#else /* sparc (and others, you figure it out) */
/* for Sun and others do something like: */
#define mGetTime(t)  times(&Timebuf) ;                               \
   (t) = Timebuf.tms_utime + Timebuf.tms_stime ;
#endif /* sparc */

#endif
#endif

#ifdef	DISPLAY
#include <starbase.c.h>		/* HP display */
#endif

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

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

#define X	0
#define Y	1

#ifdef	DISPLAY
int	Display_Tests = 0 ;
double	Display_Scale ;
double	Display_OffsetX ;
double	Display_OffsetY ;
int	Fd ;
#endif

typedef struct {
	double	time_total ;
	int	test_ratio ;
	int	test_times ;
	int	work ;
	char	*name ;
	int	flag ;
} Statistics, *pStatistics ;

#define ANGLE_TEST		 0
#define BARYCENTRIC_TEST	 1
#define CROSSINGS_TEST		 2
#define EXTERIOR_TEST		 3
#define GRID_TEST		 4
#define INCLUSION_TEST		 5
#define CROSSMULT_TEST		 6
#define PLANE_TEST		 7
#define SPACKMAN_TEST		 8
#define TRAPEZOID_TEST		 9
#define WEILER_TEST		 10
/* +++ add new name here and increment TOT_NUM_TESTS +++ */
#define TOT_NUM_TESTS		 11

Statistics	St[TOT_NUM_TESTS] ;

char	*TestName[] = {
	"angle",
	"barycentric",
	"crossings",
	"exterior",
	"grid",
	"inclusion",
	"cross-mult",
	"plane",
	"spackman",
	"trapezoid",
	"weiler" } ;
/* +++ add new name here +++ */

/* minimum & maximum number of polygon vertices to generate */
#define TOT_VERTS	1000
static int	Min_Verts = 3 ;
static int	Max_Verts = 6 ;

/* Test polygons are generated by going CCW in a circle around the origin from
 * the X+ axis around and generating vertices.	The radius is how big the
 * circumscribing circle is, the perturbation is how much each vertex is varied.
 * So, radius 1 and perturbation 0 gives a regular, inscribed polygon;
 * radius 0 and perturbation 1 gives a totally random polygon in the
 * space [-1,1)
 */
static double	Vertex_Radius = 1.0 ;
static double	Vertex_Perturbation = 0.0 ;

/* A box is circumscribed around the test polygon.  Making this box bigger
 * is useful for a higher rejection rate.  For example, a ray tracing bounding
 * box usually contains a few polygons, so making the box ratio say 2 or so
 * could simulate this type of box.
 */
static double	Box_Ratio = 1.0 ;

/* for debugging purposes, you might want to set Test_Polygons and Test_Points
 * high (say 1000), and then set the *_Test_Times to 1.	 The timings will be
 * useless, but you'll test 1000 polygons each with 1000 points.  You'll also
 * probably want to set the Vertex_Perturbation to something > 0.0.
 */
/* number of different polygons to try - I like 50 or so; left low for now */
static int	Test_Polygons = 20 ;

/* number of different intersection points to try - I like 50 or so */
static int	Test_Points = 20 ;

/* for debugging or constrained test purposes, this value constrains the value
 * of the polygon points and the test points to those on a grid emanating from
 * the origin with this increment.  Points are shifted to the closest grid
 * point.  0.0 means no grid.  NOTE:  by setting this value, a large number
 * of points will be generated exactly on interior (triangle fan) or exterior
 * edges.  Interior edge points will cause problems for the algorithms that
 * generate interior edges (triangle fan).  "On edge" points are arbitrarily
 * determined to be inside or outside the polygon, so results can differ.
 */
static double	Constraint_Increment = 0.0 ;

/* default resolutions */
static	int	Grid_Resolution = 20 ;
static	int	Trapezoid_Bins = 20 ;

#define Max(a,b)	(((a)>(b))?(a):(b))

#define FPRINTF_POLYGON						\
		fprintf( stderr, "point %g %g\n",		\
		    (float)point[X], (float)point[Y] ) ;	\
		fprintf( stderr, "polygon (%d vertices):\n",	\
			numverts ) ;				\
		for ( n = 0 ; n < numverts ; n++ ) {		\
		    fprintf( stderr, "	%g %g\n",		\
			(float)pgon[n][X], (float)pgon[n][Y]);	\
		}

/* timing functions */
#ifdef	TIMER

#define START_TIMER( test_id )						\
	/* do the test a bunch of times to get a useful time reading */ \
	mGetTime( timestart ) ;						\
	for ( tcnt = St[test_id].test_times+1 ; --tcnt ; )

#define STOP_TIMER( test_id )						\
	mGetTime( timestop ) ;						\
	/* time in microseconds */					\
	St[test_id].time_total +=					\
	    1000000.0 * (double)(timestop - timestart) /		\
	    (double)(HZ * St[test_id].test_times) ;
#else
#define START_TIMER( test_id )
#define STOP_TIMER( test_id )
#endif

char	*getenv() ;
void	Usage() ;
void	ScanOpts() ;
void	ConstrainPoint() ;
void	BreakString() ;
#ifdef	DISPLAY
void	DisplayPolygon() ;
void	DisplayPoint() ;
#endif


/* test program - see Usage() for command line options */
main(argc,argv)
int argc;  char *argv[];
{
#ifdef	TIMER
register long	tcnt ;
long	timestart ;
long	timestop ;
#endif

int	i, j, k, n, numverts, inside_flag, inside_tot, numrec ;
double	pgon[TOT_VERTS][2], point[2], angle, ran_offset ;
double	rangex, rangey, scale, minx, maxx, diffx, miny, maxy, diffy ;
double	offx, offy ;
char	str[256], *strplus ;
GridSet grid_set ;
pPlaneSet	p_plane_set ;
pSpackmanSet	p_spackman_set ;
TrapezoidSet	trap_set ;

#ifdef	CONVEX
pPlaneSet	p_ext_set ;
pInclusionAnchor	p_inc_anchor ;
#endif

    SRAN() ;

    ScanOpts( argc, argv ) ;

    for ( i = 0 ; i < TOT_NUM_TESTS ; i++ ) {
	St[i].time_total = 0.0 ;
	if ( i == ANGLE_TEST ) {
	    /* angle test is real slow, so test it fewer times */
	    St[i].test_ratio = MACHINE_TEST_RATIO / 10 ;
	} else {
	    St[i].test_ratio = MACHINE_TEST_RATIO ;
	}
	St[i].name = TestName[i] ;
	St[i].flag = 0 ;
    }

    inside_tot = 0 ;

#ifdef	CONVEX
    if ( Vertex_Perturbation > 0.0 && Max_Verts > 3 ) {
	fprintf( stderr,
		"warning: vertex perturbation is > 0.0, which is exciting\n");
	fprintf( stderr,
		"    when using convex-only algorithms!\n" ) ;
    }
#endif

    if ( Min_Verts == Max_Verts ) {
	sprintf( str, "\nPolygons with %d vertices, radius %g, "
		"perturbation +/- %g, bounding box scale %g",
	    Min_Verts, Vertex_Radius, Vertex_Perturbation, Box_Ratio ) ;
    } else {
	sprintf( str, "\nPolygons with %d to %d vertices, radius %g, "
		"perturbation +/- %g, bounding box scale %g",
	    Min_Verts, Max_Verts, Vertex_Radius, Vertex_Perturbation,
	    Box_Ratio ) ;
    }
    strplus = &str[strlen(str)] ;
    if ( St[TRAPEZOID_TEST].work ) {
	sprintf( strplus, ", %d trapezoid bins", Trapezoid_Bins ) ;
	strplus = &str[strlen(str)] ;
    }
    if ( St[GRID_TEST].work ) {
	sprintf( strplus, ", %d grid resolution", Grid_Resolution ) ;
	strplus = &str[strlen(str)] ;
    }
#ifdef	CONVEX
    sprintf( strplus, ", convex" ) ;
    strplus = &str[strlen(str)] ;
#ifdef	HYBRID
    sprintf( strplus, ", hybrid" ) ;
    strplus = &str[strlen(str)] ;
#endif
#endif
#ifdef	SORT
    if ( St[PLANE_TEST].work || St[SPACKMAN_TEST].work ) {
	sprintf( strplus, ", using triangles sorted by edge lengths" ) ;
	strplus = &str[strlen(str)] ;
#ifdef	CONVEX
	sprintf( strplus, " and areas" ) ;
	strplus = &str[strlen(str)] ;
#endif
    }
#endif
#ifdef	RANDOM
    if ( St[EXTERIOR_TEST].work ) {
	sprintf( strplus, ", exterior edges' order randomized" ) ;
	strplus = &str[strlen(str)] ;
    }
#endif
    sprintf( strplus, ".\n" ) ;
    strplus = &str[strlen(str)] ;
    BreakString( str ) ;
    printf( "%s", str ) ;

    printf(
	" Testing %d polygons with %d points\n", Test_Polygons, Test_Points ) ;

#ifdef	TIMER
    printf( "doing timings" ) ;
    fflush( stdout ) ;
#endif
    for ( i = 0 ; i < Test_Polygons ; i++ ) {

	/* make an arbitrary polygon fitting 0-1 range in x and y */
	numverts = Min_Verts +
		(int)(RAN01() * (double)(Max_Verts-Min_Verts+1)) ;

	/* add a random offset to the angle so that each polygon is not in
	 * some favorable (or unfavorable) alignment.
	 */
	ran_offset = 2.0 * M_PI * RAN01() ;
	minx = miny =  99999.0 ;
	maxx = maxy = -99999.0 ;
	for ( j = 0 ; j < numverts ; j++ ) {
	    angle = 2.0 * M_PI * (double)j / (double)numverts + ran_offset ;
	    pgon[j][X] = cos(angle) * Vertex_Radius +
		( RAN01() * 2.0 - 1.0 ) * Vertex_Perturbation ;
	    pgon[j][Y] = sin(angle) * Vertex_Radius +
		( RAN01() * 2.0 - 1.0 ) * Vertex_Perturbation ;

	    ConstrainPoint( pgon[j] ) ;

	    if ( pgon[j][X] < minx ) minx = pgon[j][X] ;
	    if ( pgon[j][X] > maxx ) maxx = pgon[j][X] ;
	    if ( pgon[j][Y] < miny ) miny = pgon[j][Y] ;
	    if ( pgon[j][Y] > maxy ) maxy = pgon[j][Y] ;
	}

	offx = ( maxx + minx ) / 2.0 ;
	offy = ( maxy + miny ) / 2.0 ;
	if ( (diffx = maxx - minx ) > ( diffy = maxy - miny ) ) {
	    scale = 2.0 / (Box_Ratio * diffx) ;
	    rangex = 1.0 ;
	    rangey = diffy / diffx ;
	} else {
	    scale = 2.0 / (Box_Ratio * diffy) ;
	    rangex = diffx / diffy ;
	    rangey = 1.0 ;
	}

	for ( j = 0 ; j < numverts ; j++ ) {
	    pgon[j][X] = ( pgon[j][X] - offx ) * scale ;
	    pgon[j][Y] = ( pgon[j][Y] - offy ) * scale ;
	}

	/* Set up number of times to test a point against a polygon, for
	 * the sake of getting a reasonable timing.  We already know how
	 * most of these will perform, so scale their # tests accordingly.
	 */
	for ( j = 0 ; j < TOT_NUM_TESTS ; j++ ) {
	    if ( ( j == GRID_TEST ) || ( j == TRAPEZOID_TEST ) ) {
		St[j].test_times = Max( St[j].test_ratio /
		    (int)sqrt((double)numverts), 1 ) ;
	    } else {
		St[j].test_times = Max( St[j].test_ratio / numverts, 1 ) ;
	    }
	}

	/* set up tests */
#ifdef	CONVEX
	if ( St[EXTERIOR_TEST].work ) {
	    p_ext_set = ExteriorSetup( pgon, numverts ) ;
	}
#endif

	if ( St[GRID_TEST].work ) {
	    GridSetup( pgon, numverts, Grid_Resolution, &grid_set ) ;
	}

#ifdef	CONVEX
	if ( St[INCLUSION_TEST].work ) {
	    p_inc_anchor = InclusionSetup( pgon, numverts ) ;
	}
#endif

	if ( St[PLANE_TEST].work ) {
	    p_plane_set = PlaneSetup( pgon, numverts ) ;
	}

	if ( St[SPACKMAN_TEST].work ) {
	    p_spackman_set = SpackmanSetup( pgon, numverts, &numrec ) ;
	}

	if ( St[TRAPEZOID_TEST].work ) {
	    TrapezoidSetup( pgon, numverts, Trapezoid_Bins, &trap_set ) ;
	}

#ifdef	DISPLAY
	if ( Display_Tests ) {
	    DisplayPolygon( pgon, numverts, i ) ;
	}
#endif

	/* now try # of points against it */
	for ( j = 0 ; j < Test_Points ; j++ ) {
	    point[X] = RAN01() * rangex * 2.0 - rangex ;
	    point[Y] = RAN01() * rangey * 2.0 - rangey ;

	    ConstrainPoint( point ) ;

#ifdef	DISPLAY
	    if ( Display_Tests ) {
		DisplayPoint( point, TRUE ) ;
	    }
#endif

	    if ( St[ANGLE_TEST].work ) {
		START_TIMER( ANGLE_TEST )
		    St[ANGLE_TEST].flag = AngleTest( pgon, numverts, point ) ;
		STOP_TIMER( ANGLE_TEST )
	    }
	    if ( St[BARYCENTRIC_TEST].work ) {
		START_TIMER( BARYCENTRIC_TEST )
		    St[BARYCENTRIC_TEST].flag =
			    BarycentricTest( pgon, numverts, point ) ;
		STOP_TIMER( BARYCENTRIC_TEST )
	    }
	    if ( St[CROSSINGS_TEST].work ) {
		START_TIMER( CROSSINGS_TEST )
		    St[CROSSINGS_TEST].flag =
			    CrossingsTest( pgon, numverts, point ) ;
		STOP_TIMER( CROSSINGS_TEST )
	    }
#ifdef	CONVEX
	    if ( St[EXTERIOR_TEST].work ) {
		START_TIMER( EXTERIOR_TEST )
		    St[EXTERIOR_TEST].flag =
			    ExteriorTest( p_ext_set, numverts, point );
		STOP_TIMER( EXTERIOR_TEST )
	    }
#endif
	    if ( St[GRID_TEST].work ) {
		START_TIMER( GRID_TEST )
		    St[GRID_TEST].flag = GridTest( &grid_set, point ) ;
		STOP_TIMER( GRID_TEST )
	    }
#ifdef	CONVEX
	    if ( St[INCLUSION_TEST].work ) {
		START_TIMER( INCLUSION_TEST )
		    St[INCLUSION_TEST].flag =
			    InclusionTest( p_inc_anchor, point );
		STOP_TIMER( INCLUSION_TEST )
	    }
#endif
	    if ( St[CROSSMULT_TEST].work ) {
		START_TIMER( CROSSMULT_TEST )
		    St[CROSSMULT_TEST].flag = CrossingsMultiplyTest(
				pgon, numverts, point ) ;