Graphics Gems 1~5全系列的代码分享

/* ptinpoly.c - point in polygon inside/outside code. by Eric Haines, 3D/Eye Inc, erich@eye.com 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 ===