home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsiv / ptpoly_haines / ptinpoly.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-12-13  |  55.9 KB  |  2,139 lines
  1. /* ptinpoly.c - point in polygon inside/outside code.
  2.  
  3.    by Eric Haines, 3D/Eye Inc, erich@eye.com
  4.  
  5.    This code contains the following algorithms:
  6.     crossings - count the crossing made by a ray from the test point
  7.     crossings-multiply - as above, but avoids a division; often a bit faster
  8.     angle summation - sum the angle formed by point and vertex pairs
  9.     weiler angle summation - sum the angles using quad movements
  10.     half-plane testing - test triangle fan using half-space planes
  11.     barycentric coordinates - test triangle fan w/barycentric coords
  12.     spackman barycentric - preprocessed barycentric coordinates
  13.     trapezoid testing - bin sorting algorithm
  14.     grid testing - grid imposed on polygon
  15.     exterior test - for convex polygons, check exterior of polygon
  16.     inclusion test - for convex polygons, use binary search for edge.
  17. */
  18.  
  19. #include <stdio.h>
  20. #include <stdlib.h>
  21. #include <math.h>
  22. #include "ptinpoly.h"
  23.  
  24. #define X    0
  25. #define Y    1
  26.  
  27. #ifndef TRUE
  28. #define TRUE    1
  29. #define FALSE    0
  30. #endif
  31.  
  32. #ifndef HUGE
  33. #define HUGE    1.797693134862315e+308
  34. #endif
  35.  
  36. #ifndef M_PI
  37. #define M_PI    3.14159265358979323846
  38. #endif
  39.  
  40. /* test if a & b are within epsilon.  Favors cases where a < b */
  41. #define Near(a,b,eps)    ( ((b)-(eps)<(a)) && ((a)-(eps)<(b)) )
  42.  
  43. #define MALLOC_CHECK( a )    if ( !(a) ) {                   \
  44.                     fprintf( stderr, "out of memory\n" ) ; \
  45.                     exit(1) ;                   \
  46.                 }
  47.  
  48.  
  49. /* ======= Crossings algorithm ============================================ */
  50.  
  51. /* Shoot a test ray along +X axis.  The strategy, from MacMartin, is to
  52.  * compare vertex Y values to the testing point's Y and quickly discard
  53.  * edges which are entirely to one side of the test ray.
  54.  *
  55.  * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
  56.  * _point_, returns 1 if inside, 0 if outside.    WINDING and CONVEX can be
  57.  * defined for this test.
  58.  */
  59. int CrossingsTest( pgon, numverts, point )
  60. double    pgon[][2] ;
  61. int    numverts ;
  62. double    point[2] ;
  63. {
  64. #ifdef    WINDING
  65. register int    crossings ;
  66. #endif
  67. register int    j, yflag0, yflag1, inside_flag, xflag0 ;
  68. register double ty, tx, *vtx0, *vtx1 ;
  69. #ifdef    CONVEX
  70. register int    line_flag ;
  71. #endif
  72.  
  73.     tx = point[X] ;
  74.     ty = point[Y] ;
  75.  
  76.     vtx0 = pgon[numverts-1] ;
  77.     /* get test bit for above/below X axis */
  78.     yflag0 = ( vtx0[Y] >= ty ) ;
  79.     vtx1 = pgon[0] ;
  80.  
  81. #ifdef    WINDING
  82.     crossings = 0 ;
  83. #else
  84.     inside_flag = 0 ;
  85. #endif
  86. #ifdef    CONVEX
  87.     line_flag = 0 ;
  88. #endif
  89.     for ( j = numverts+1 ; --j ; ) {
  90.  
  91.     yflag1 = ( vtx1[Y] >= ty ) ;
  92.     /* check if endpoints straddle (are on opposite sides) of X axis
  93.      * (i.e. the Y's differ); if so, +X ray could intersect this edge.
  94.      */
  95.     if ( yflag0 != yflag1 ) {
  96.         xflag0 = ( vtx0[X] >= tx ) ;
  97.         /* check if endpoints are on same side of the Y axis (i.e. X's
  98.          * are the same); if so, it's easy to test if edge hits or misses.
  99.          */
  100.         if ( xflag0 == ( vtx1[X] >= tx ) ) {
  101.  
  102.         /* if edge's X values both right of the point, must hit */
  103. #ifdef    WINDING
  104.         if ( xflag0 ) crossings += ( yflag0 ? -1 : 1 ) ;
  105. #else
  106.         if ( xflag0 ) inside_flag = !inside_flag ;
  107. #endif
  108.         } else {
  109.         /* compute intersection of pgon segment with +X ray, note
  110.          * if >= point's X; if so, the ray hits it.
  111.          */
  112.         if ( (vtx1[X] - (vtx1[Y]-ty)*
  113.              ( vtx0[X]-vtx1[X])/(vtx0[Y]-vtx1[Y])) >= tx ) {
  114. #ifdef    WINDING
  115.             crossings += ( yflag0 ? -1 : 1 ) ;
  116. #else
  117.             inside_flag = !inside_flag ;
  118. #endif
  119.         }
  120.         }
  121. #ifdef    CONVEX
  122.         /* if this is second edge hit, then done testing */
  123.         if ( line_flag ) goto Exit ;
  124.  
  125.         /* note that one edge has been hit by the ray's line */
  126.         line_flag = TRUE ;
  127. #endif
  128.     }
  129.  
  130.     /* move to next pair of vertices, retaining info as possible */
  131.     yflag0 = yflag1 ;
  132.     vtx0 = vtx1 ;
  133.     vtx1 += 2 ;
  134.     }
  135. #ifdef    CONVEX
  136.     Exit: ;
  137. #endif
  138. #ifdef    WINDING
  139.     /* test if crossings is not zero */
  140.     inside_flag = (crossings != 0) ;
  141. #endif
  142.  
  143.     return( inside_flag ) ;
  144. }
  145.  
  146. /* ======= Angle summation algorithm ======================================= */
  147.  
  148. /* Sum angles made by (vtxN to test point to vtxN+1), check if in proper
  149.  * range to be inside.    VERY SLOW, included for tutorial reasons (i.e.
  150.  * to show why one should never use this algorithm).
  151.  *
  152.  * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
  153.  * _point_, returns 1 if inside, 0 if outside.
  154.  */
  155. int AngleTest( pgon, numverts, point )
  156. double    pgon[][2] ;
  157. int    numverts ;
  158. double    point[2] ;
  159. {
  160. register double *vtx0, *vtx1, angle, len, vec0[2], vec1[2], vec_dot ;
  161. register int    j ;
  162. int    inside_flag ;
  163.  
  164.     /* sum the angles and see if answer mod 2*PI > PI */
  165.     vtx0 = pgon[numverts-1] ;
  166.     vec0[X] = vtx0[X] - point[X] ;
  167.     vec0[Y] = vtx0[Y] - point[Y] ;
  168.     len = sqrt( vec0[X] * vec0[X] + vec0[Y] * vec0[Y] ) ;
  169.     if ( len <= 0.0 ) {
  170.     /* point and vertex coincide */
  171.     return( 1 ) ;
  172.     }
  173.     vec0[X] /= len ;
  174.     vec0[Y] /= len ;
  175.  
  176.     angle = 0.0 ;
  177.     for ( j = 0 ; j < numverts ; j++ ) {
  178.     vtx1 = pgon[j] ;
  179.     vec1[X] = vtx1[X] - point[X] ;
  180.     vec1[Y] = vtx1[Y] - point[Y] ;
  181.     len = sqrt( vec1[X] * vec1[X] + vec1[Y] * vec1[Y] ) ;
  182.     if ( len <= 0.0 ) {
  183.         /* point and vertex coincide */
  184.         return( 1 ) ;
  185.     }
  186.     vec1[X] /= len ;
  187.     vec1[Y] /= len ;
  188.  
  189.     /* check if vec1 is to "left" or "right" of vec0 */
  190.     vec_dot = vec0[X] * vec1[X] + vec0[Y] * vec1[Y] ;
  191.     if ( vec_dot < -1.0 ) {
  192.         /* point is on edge, so always add 180 degrees.  always
  193.          * adding is not necessarily the right answer for
  194.          * concave polygons and subtractive triangles.
  195.          */
  196.         angle += M_PI ;
  197.     } else if ( vec_dot < 1.0 ) {
  198.         if ( vec0[X] * vec1[Y] - vec1[X] * vec0[Y] >= 0.0 ) {
  199.         /* add angle due to dot product of vectors */
  200.         angle += acos( vec_dot ) ;
  201.         } else {
  202.         /* subtract angle due to dot product of vectors */
  203.         angle -= acos( vec_dot ) ;
  204.         }
  205.     } /* if vec_dot >= 1.0, angle does not change */
  206.  
  207.     /* get to next point */
  208.     vtx0 = vtx1 ;
  209.     vec0[X] = vec1[X] ;
  210.     vec0[Y] = vec1[Y] ;
  211.     }
  212.     /* test if between PI and 3*PI, 5*PI and 7*PI, etc */
  213.     inside_flag = fmod( fabs(angle) + M_PI, 4.0*M_PI ) > 2.0*M_PI ;
  214.  
  215.     return( inside_flag ) ;
  216. }
  217.  
  218. /* ======= Weiler algorithm ============================================ */
  219.  
  220. /* Track quadrant movements using Weiler's algorithm (elsewhere in Graphics
  221.  * Gems IV).  Algorithm has been optimized for testing purposes, but the
  222.  * crossings algorithm is still faster.     Included to provide timings.
  223.  *
  224.  * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
  225.  * _point_, returns 1 if inside, 0 if outside.    WINDING can be defined for
  226.  * this test.
  227.  */
  228.  
  229. #define QUADRANT( vtx, x, y )    \
  230.     (((vtx)[X]>(x)) ? ( ((vtx)[Y]>(y)) ? 0:3 ) : ( ((vtx)[Y]>(y)) ? 1:2 ))
  231.  
  232. #define X_INTERCEPT( v0, v1, y )    \
  233.     ( (((v1)[X]-(v0)[X])/((v1)[Y]-(v0)[Y])) * ((y)-(v0)[Y]) + (v0)[X] )
  234.  
  235. int WeilerTest( pgon, numverts, point )
  236. double    pgon[][2] ;
  237. int    numverts ;
  238. double    point[2] ;
  239. {
  240. register int    angle, qd, next_qd, delta, j ;
  241. register double ty, tx, *vtx0, *vtx1 ;
  242. int    inside_flag ;
  243.  
  244.     tx = point[X] ;
  245.     ty = point[Y] ;
  246.  
  247.     vtx0 = pgon[numverts-1] ;
  248.     qd = QUADRANT( vtx0, tx, ty ) ;
  249.     angle = 0 ;
  250.  
  251.     vtx1 = pgon[0] ;
  252.  
  253.     for ( j = numverts+1 ; --j ; ) {
  254.     /* calculate quadrant and delta from last quadrant */
  255.     next_qd = QUADRANT( vtx1, tx, ty ) ;
  256.     delta = next_qd - qd ;
  257.  
  258.     /* adjust delta and add it to total angle sum */
  259.     switch ( delta ) {
  260.         case 0:    /* do nothing */
  261.         break ;
  262.         case -1:
  263.         case 3:
  264.         angle-- ;
  265.         qd = next_qd ;
  266.         break ;
  267.  
  268.         case 1:
  269.         case -3:
  270.         angle++ ;
  271.         qd = next_qd ;
  272.         break ;
  273.  
  274.         case 2:
  275.         case -2:
  276.         if (X_INTERCEPT( vtx0, vtx1, ty ) > tx ) {
  277.             angle -= delta ;
  278.         } else {
  279.             angle += delta ;
  280.         }
  281.         qd = next_qd ;
  282.         break ;
  283.     }
  284.  
  285.     /* increment for next step */
  286.     vtx0 = vtx1 ;
  287.     vtx1 += 2 ;
  288.     }
  289.  
  290. #ifdef    WINDING
  291.     /* simple windings test:  if angle != 0, then point is inside */
  292.     inside_flag = ( angle != 0 ) ;
  293. #else
  294.     /* Jordan test:  if angle is +-4, 12, 20, etc then point is inside */
  295.     inside_flag = ( (angle/4) & 0x1 ) ;
  296. #endif
  297.  
  298.     return (inside_flag) ;
  299. }
  300.  
  301. #undef    QUADRANT
  302. #undef    Y_INTERCEPT
  303.  
  304. /* ======= Triangle half-plane algorithm ================================== */
  305.  
  306. /* Split the polygon into a fan of triangles and for each triangle test if
  307.  * the point is inside of the three half-planes formed by the triangle's edges.
  308.  *
  309.  * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
  310.  * which returns a pointer to a plane set array.
  311.  * Call testing procedure with a pointer to this array, _numverts_, and
  312.  * test point _point_, returns 1 if inside, 0 if outside.
  313.  * Call cleanup with pointer to plane set array to free space.
  314.  *
  315.  * SORT and CONVEX can be defined for this test.
  316.  */
  317.  
  318. /* split polygons along set of x axes - call preprocess once */
  319. pPlaneSet    PlaneSetup( pgon, numverts )
  320. double    pgon[][2] ;
  321. int    numverts ;
  322. {
  323. int    i, p1, p2 ;
  324. double    tx, ty, vx0, vy0 ;
  325. pPlaneSet    pps, pps_return ;
  326. #ifdef    SORT
  327. double    len[3], len_temp ;
  328. int    j ;
  329. PlaneSet    ps_temp ;
  330. #ifdef    CONVEX
  331. pPlaneSet    pps_new ;
  332. pSizePlanePair p_size_pair ;
  333. #endif
  334. #endif
  335.  
  336.     pps = pps_return =
  337.         (pPlaneSet)malloc( 3 * (numverts-2) * sizeof( PlaneSet )) ;
  338.     MALLOC_CHECK( pps ) ;
  339. #ifdef    CONVEX
  340. #ifdef    SORT
  341.     p_size_pair =
  342.     (pSizePlanePair)malloc( (numverts-2) * sizeof( SizePlanePair ) ) ;
  343.     MALLOC_CHECK( p_size_pair ) ;
  344. #endif
  345. #endif
  346.  
  347.     vx0 = pgon[0][X] ;
  348.     vy0 = pgon[0][Y] ;
  349.  
  350.     for ( p1 = 1, p2 = 2 ; p2 < numverts ; p1++, p2++ ) {
  351.     pps->vx = vy0 - pgon[p1][Y] ;
  352.     pps->vy = pgon[p1][X] - vx0 ;
  353.     pps->c = pps->vx * vx0 + pps->vy * vy0 ;
  354. #ifdef    SORT
  355.     len[0] = pps->vx * pps->vx + pps->vy * pps->vy ;
  356. #ifdef    CONVEX
  357. #ifdef    HYBRID
  358.     pps->ext_flag = ( p1 == 1 ) ;
  359. #endif
  360.     /* Sort triangles by areas, so compute (twice) the area here */
  361.     p_size_pair[p1-1].pps = pps ;
  362.     p_size_pair[p1-1].size =
  363.             ( pgon[0][X] * pgon[p1][Y] ) +
  364.             ( pgon[p1][X] * pgon[p2][Y] ) +
  365.             ( pgon[p2][X] * pgon[0][Y] ) -
  366.             ( pgon[p1][X] * pgon[0][Y] ) -
  367.             ( pgon[p2][X] * pgon[p1][Y] ) -
  368.             ( pgon[0][X] * pgon[p2][Y] ) ;
  369. #endif
  370. #endif
  371.     pps++ ;
  372.     pps->vx = pgon[p1][Y] - pgon[p2][Y] ;
  373.     pps->vy = pgon[p2][X] - pgon[p1][X] ;
  374.     pps->c = pps->vx * pgon[p1][X] + pps->vy * pgon[p1][Y] ;
  375. #ifdef    SORT
  376.     len[1] = pps->vx * pps->vx + pps->vy * pps->vy ;
  377. #endif
  378. #ifdef    CONVEX
  379. #ifdef    HYBRID
  380.     pps->ext_flag = TRUE ;
  381. #endif
  382. #endif
  383.     pps++ ;
  384.     pps->vx = pgon[p2][Y] - vy0 ;
  385.     pps->vy = vx0 - pgon[p2][X] ;
  386.     pps->c = pps->vx * pgon[p2][X] + pps->vy * pgon[p2][Y] ;
  387. #ifdef    SORT
  388.     len[2] = pps->vx * pps->vx + pps->vy * pps->vy ;
  389. #endif
  390. #ifdef    CONVEX
  391. #ifdef    HYBRID
  392.     pps->ext_flag = ( p2 == numverts-1 ) ;
  393. #endif
  394. #endif
  395.  
  396.     /* find an average point which must be inside of the triangle */
  397.     tx = ( vx0 + pgon[p1][X] + pgon[p2][X] ) / 3.0 ;
  398.     ty = ( vy0 + pgon[p1][Y] + pgon[p2][Y] ) / 3.0 ;
  399.  
  400.     /* check sense and reverse if test point is not thought to be inside
  401.      * first triangle
  402.      */
  403.     if ( pps->vx * tx + pps->vy * ty >= pps->c ) {
  404.         /* back up to start of plane set */
  405.         pps -= 2 ;
  406.         /* point is thought to be outside, so reverse sense of edge
  407.          * normals so that it is correctly considered inside.
  408.          */
  409.         for ( i = 0 ; i < 3 ; i++ ) {
  410.         pps->vx = -pps->vx ;
  411.         pps->vy = -pps->vy ;
  412.         pps->c    = -pps->c ;
  413.         pps++ ;
  414.         }
  415.     } else {
  416.         pps++ ;
  417.     }
  418.  
  419. #ifdef    SORT
  420.     /* sort the planes based on the edge lengths */
  421.     pps -= 3 ;
  422.     for ( i = 0 ; i < 2 ; i++ ) {
  423.         for ( j = i+1 ; j < 3 ; j++ ) {
  424.         if ( len[i] < len[j] ) {
  425.             ps_temp = pps[i] ;
  426.             pps[i] = pps[j] ;
  427.             pps[j] = ps_temp ;
  428.             len_temp = len[i] ;
  429.             len[i] = len[j] ;
  430.             len[j] = len_temp ;
  431.         }
  432.         }
  433.     }
  434.     pps += 3 ;
  435. #endif
  436.     }
  437.  
  438. #ifdef    CONVEX
  439. #ifdef    SORT
  440.     /* sort the triangles based on their areas */
  441.     qsort( p_size_pair, numverts-2,
  442.         sizeof( SizePlanePair ), CompareSizePlanePairs ) ;
  443.  
  444.     /* make the plane sets match the sorted order */
  445.     for ( i = 0, pps = pps_return
  446.     ; i < numverts-2
  447.     ; i++ ) {
  448.  
  449.     pps_new = p_size_pair[i].pps ;
  450.     for ( j = 0 ; j < 3 ; j++, pps++, pps_new++ ) {
  451.         ps_temp = *pps ;
  452.         *pps = *pps_new ;
  453.         *pps_new = ps_temp ;
  454.     }
  455.     }
  456.     free( p_size_pair ) ;
  457. #endif
  458. #endif
  459.  
  460.     return( pps_return ) ;
  461. }
  462.  
  463. #ifdef    CONVEX
  464. #ifdef    SORT
  465. int CompareSizePlanePairs( p_sp0, p_sp1 )
  466. pSizePlanePair    p_sp0, p_sp1 ;
  467. {
  468.     if ( p_sp0->size == p_sp1->size ) {
  469.     return( 0 ) ;
  470.     } else {
  471.     return( p_sp0->size > p_sp1->size ? -1 : 1 ) ;
  472.     }
  473. }
  474. #endif
  475. #endif
  476.  
  477.  
  478. /* check point for inside of three "planes" formed by triangle edges */
  479. int PlaneTest( p_plane_set, numverts, point )
  480. pPlaneSet    p_plane_set ;
  481. int    numverts ;
  482. double    point[2] ;
  483. {
  484. register pPlaneSet    ps ;
  485. register int    p2 ;
  486. #ifndef CONVEX
  487. register int    inside_flag ;
  488. #endif
  489. register double tx, ty ;
  490.  
  491.     tx = point[X] ;
  492.     ty = point[Y] ;
  493.  
  494. #ifndef CONVEX
  495.     inside_flag = 0 ;
  496. #endif
  497.  
  498.     for ( ps = p_plane_set, p2 = numverts-1 ; --p2 ; ) {
  499.  
  500.     if ( ps->vx * tx + ps->vy * ty < ps->c ) {
  501.         ps++ ;
  502.         if ( ps->vx * tx + ps->vy * ty < ps->c ) {
  503.         ps++ ;
  504.         /* note: we make the third edge have a slightly different
  505.          * equality condition, since this third edge is in fact
  506.          * the next triangle's first edge.  Not fool-proof, but
  507.          * it doesn't hurt (better would be to keep track of the
  508.          * triangle's area sign so we would know which kind of
  509.          * triangle this is).  Note that edge sorting nullifies
  510.          * this special inequality, too.
  511.          */
  512.         if ( ps->vx * tx + ps->vy * ty <= ps->c ) {
  513.             /* point is inside polygon */
  514. #ifdef CONVEX
  515.             return( 1 ) ;
  516. #else
  517.             inside_flag = !inside_flag ;
  518. #endif
  519.         }
  520. #ifdef    CONVEX
  521. #ifdef    HYBRID
  522.         /* check if outside exterior edge */
  523.         else if ( ps->ext_flag ) return( 0 ) ;
  524. #endif
  525. #endif
  526.         ps++ ;
  527.         } else {
  528. #ifdef    CONVEX
  529. #ifdef    HYBRID
  530.         /* check if outside exterior edge */
  531.         if ( ps->ext_flag ) return( 0 ) ;
  532. #endif
  533. #endif
  534.         /* get past last two plane tests */
  535.         ps += 2 ;
  536.         }
  537.     } else {
  538. #ifdef    CONVEX
  539. #ifdef    HYBRID
  540.         /* check if outside exterior edge */
  541.         if ( ps->ext_flag ) return( 0 ) ;
  542. #endif
  543. #endif
  544.         /* get past all three plane tests */
  545.         ps += 3 ;
  546.     }
  547.     }
  548.  
  549. #ifdef CONVEX
  550.     /* for convex, if we make it to here, all triangles were missed */
  551.     return( 0 ) ;
  552. #else
  553.     return( inside_flag ) ;
  554. #endif
  555. }
  556.  
  557. void PlaneCleanup( p_plane_set )
  558. pPlaneSet    p_plane_set ;
  559. {
  560.     free( p_plane_set ) ;
  561. }
  562.  
  563. /* ======= Barycentric algorithm ========================================== */
  564.  
  565. /* Split the polygon into a fan of triangles and for each triangle test if
  566.  * the point has barycentric coordinates which are inside the triangle.
  567.  * Similar to Badouel's code in Graphics Gems, with a little more efficient
  568.  * coding.
  569.  *
  570.  * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
  571.  * _point_, returns 1 if inside, 0 if outside.
  572.  */
  573. int BarycentricTest( pgon, numverts, point )
  574. double    pgon[][2] ;
  575. int    numverts ;
  576. double    point[2] ;
  577. {
  578. register double *pg1, *pg2, *pgend ;
  579. register double tx, ty, u0, u1, u2, v0, v1, vx0, vy0, alpha, beta, denom ;
  580. int    inside_flag ;
  581.  
  582.     tx = point[X] ;
  583.     ty = point[Y] ;
  584.     vx0 = pgon[0][X] ;
  585.     vy0 = pgon[0][Y] ;
  586.     u0 = tx - vx0 ;
  587.     v0 = ty - vy0 ;
  588.  
  589.     inside_flag = 0 ;
  590.     pgend = pgon[numverts-1] ;
  591.     for ( pg1 = pgon[1], pg2 = pgon[2] ; pg1 != pgend ; pg1+=2, pg2+=2 ) {
  592.  
  593.     u1 = pg1[X] - vx0 ;
  594.     if ( u1 == 0.0 ) {
  595.  
  596.         /* 0 and 1 vertices have same X value */
  597.  
  598.         /* zero area test - can be removed for convex testing */
  599.         u2 = pg2[X] - vx0 ;
  600.         if ( ( u2 == 0.0 ) ||
  601.  
  602.         /* compute beta and check bounds */
  603.         /* we use "<= 0.0" so that points on the shared interior
  604.          * edge will (generally) be inside only one polygon.
  605.          */
  606.          ( ( beta = u0 / u2 ) <= 0.0 ) ||
  607.          ( beta > 1.0 ) ||
  608.  
  609.          /* zero area test - remove for convex testing */
  610.          ( ( v1 = pg1[Y] - vy0 ) == 0.0 ) ||
  611.  
  612.          /* compute alpha and check bounds */
  613.          ( ( alpha = ( v0 - beta *
  614.             ( pg2[Y] - vy0 ) ) / v1 ) < 0.0 ) ) {
  615.  
  616.         /* whew! missed! */
  617.         goto NextTri ;
  618.         }
  619.  
  620.     } else {
  621.         /* 0 and 1 vertices have different X value */
  622.  
  623.         /* compute denom and check for zero area triangle - check
  624.          * is not needed for convex polygon testing
  625.          */
  626.         u2 = pg2[X] - vx0 ;
  627.         v1 = pg1[Y] - vy0 ;
  628.         denom = ( pg2[Y] - vy0 ) * u1 - u2 * v1 ;
  629.         if ( ( denom == 0.0 ) ||
  630.  
  631.         /* compute beta and check bounds */
  632.         /* we use "<= 0.0" so that points on the shared interior
  633.          * edge will (generally) be inside only one polygon.
  634.          */
  635.          ( ( beta = ( v0 * u1 - u0 * v1 ) / denom ) <= 0.0 ) ||
  636.          ( beta > 1.0 ) ||
  637.  
  638.          /* compute alpha & check bounds */
  639.          ( ( alpha = ( u0 - beta * u2 ) / u1 ) < 0.0 ) ) {
  640.  
  641.         /* whew! missed! */
  642.         goto NextTri ;
  643.         }
  644.     }
  645.  
  646.     /* check gamma */
  647.     if ( alpha + beta <= 1.0 ) {
  648.         /* survived */
  649.         inside_flag = !inside_flag ;
  650.     }
  651.  
  652.     NextTri: ;
  653.     }
  654.     return( inside_flag ) ;
  655. }
  656.  
  657. /* ======= Barycentric precompute (Spackman) algorithm ===================== */
  658.  
  659. /* Split the polygon into a fan of triangles and for each triangle test if
  660.  * the point has barycentric coordinates which are inside the triangle.
  661.  * Use Spackman's normalization method to precompute various parameters.
  662.  *
  663.  * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
  664.  * which returns a pointer to the array of the parameters records and the
  665.  * number of parameter records created.
  666.  * Call testing procedure with the first vertex in the polygon _pgon[0]_,
  667.  * a pointer to this array, the number of parameter records, and test point
  668.  * _point_, returns 1 if inside, 0 if outside.
  669.  * Call cleanup with pointer to parameter record array to free space.
  670.  *
  671.  * SORT can be defined for this test.
  672.  * (CONVEX could be added: see PlaneSetup and PlaneTest for method)
  673.  */
  674. pSpackmanSet    SpackmanSetup( pgon, numverts, p_numrec )
  675. double    pgon[][2] ;
  676. int    numverts ;
  677. int    *p_numrec ;
  678. {
  679. int    p1, p2, degen ;
  680. double    denom, u1, v1, *pv[3] ;
  681. pSpackmanSet    pss, pss_return ;
  682. #ifdef    SORT
  683. double    u[2], v[2], len[2], *pv_temp ;
  684. #endif
  685.  
  686.     pss = pss_return =
  687.         (pSpackmanSet)malloc( (numverts-2) * sizeof( SpackmanSet )) ;
  688.     MALLOC_CHECK( pss ) ;
  689.  
  690.     degen = 0 ;
  691.  
  692.     for ( p1 = 1, p2 = 2 ; p2 < numverts ; p1++, p2++ ) {
  693.  
  694.     pv[0] = pgon[0] ;
  695.     pv[1] = pgon[p1] ;
  696.     pv[2] = pgon[p2] ;
  697.  
  698. #ifdef    SORT
  699.     /* Note that sorting can cause a mismatch of alpha/beta inequality
  700.      * tests.  In other words, test points on an interior line between
  701.      * test triangles will often then be wrong.
  702.      */
  703.     u[0] = pv[1][X] - pv[0][X] ;
  704.     u[1] = pv[2][X] - pv[0][X] ;
  705.     v[0] = pv[1][Y] - pv[0][Y] ;
  706.     v[1] = pv[2][Y] - pv[0][Y] ;
  707.     len[0] = u[0] * u[0] + v[0] * v[0] ;
  708.     len[1] = u[1] * u[1] + v[1] * v[1] ;
  709.  
  710.     /* compare two edges touching anchor point and put longest first */
  711.     /* we don't sort all three edges because the anchor point and
  712.      * values computed from it gets used for all triangles in the fan.
  713.      */
  714.     if ( len[0] < len[1] ) {
  715.         pv_temp = pv[1] ;
  716.         pv[1] = pv[2] ;
  717.         pv[2] = pv_temp ;
  718.     }
  719. #endif
  720.  
  721.     u1 = pv[1][X] - pv[0][X] ;
  722.     pss->u2 = pv[2][X] - pv[0][X] ;
  723.     v1 = pv[1][Y] - pv[0][Y] ;
  724.     pss->v2 = pv[2][Y] - pv[0][Y] ;
  725.     pss->u1_nonzero = !( u1 == 0.0 ) ;
  726.     if ( pss->u1_nonzero ) {
  727.         /* not zero, so compute inverse */
  728.         pss->inv_u1 = 1.0 / u1 ;
  729.         denom = pss->v2 * u1 - pss->u2 * v1 ;
  730.         if ( denom == 0.0 ) {
  731.         /* degenerate triangle, ignore it */
  732.         degen++ ;
  733.         goto Skip ;
  734.         } else {
  735.         pss->u1p = u1 / denom ;
  736.         pss->v1p = v1 / denom ;
  737.         }
  738.     } else {
  739.         if ( pss->u2 == 0.0 ) {
  740.         /* degenerate triangle, ignore it */
  741.         degen++ ;
  742.         goto Skip ;
  743.         } else {
  744.         /* not zero, so compute inverse */
  745.         pss->inv_u2 = 1.0 / pss->u2 ;
  746.         if ( v1 == 0.0 ) {
  747.             /* degenerate triangle, ignore it */
  748.             degen++ ;
  749.             goto Skip ;
  750.         } else {
  751.             pss->inv_v1 = 1.0 / v1 ;
  752.         }
  753.         }
  754.     }
  755.  
  756.     pss++ ;
  757.     Skip: ;
  758.     }
  759.  
  760.     /* number of Spackman records */
  761.     *p_numrec = numverts - degen - 2 ;
  762.     if ( degen ) {
  763.     pss = pss_return =
  764.         (pSpackmanSet)realloc( pss_return,
  765.             (numverts-2-degen) * sizeof( SpackmanSet )) ;
  766.     }
  767.  
  768.     return( pss_return ) ;
  769. }
  770.  
  771. /* barycentric, a la Gems I and Spackman's normalization precompute */
  772. int SpackmanTest( anchor, p_spackman_set, numrec, point )
  773. double    anchor[2] ;
  774. pSpackmanSet    p_spackman_set ;
  775. int    numrec ;
  776. double    point[2] ;
  777. {
  778. register pSpackmanSet    pss ;
  779. register int    inside_flag ;
  780. register int    nr ;
  781. register double tx, ty, vx0, vy0, u0, v0, alpha, beta ;
  782.  
  783.     tx = point[X] ;
  784.     ty = point[Y] ;
  785.     /* note that we really need only the first vertex of the polygon,
  786.      * so do not really need to keep the whole polygon around.
  787.      */
  788.     vx0 = anchor[X] ;
  789.     vy0 = anchor[Y] ;
  790.     u0 = tx - vx0 ;
  791.     v0 = ty - vy0 ;
  792.  
  793.     inside_flag = 0 ;
  794.  
  795.     for ( pss = p_spackman_set, nr = numrec+1 ; --nr ; pss++ ) {
  796.  
  797.     if ( pss->u1_nonzero ) {
  798.         /* 0 and 2 vertices have different X value */
  799.  
  800.         /* compute beta and check bounds */
  801.         /* we use "<= 0.0" so that points on the shared interior edge
  802.          * will (generally) be inside only one polygon.
  803.          */
  804.         beta = ( v0 * pss->u1p - u0 * pss->v1p ) ;
  805.         if ( ( beta <= 0.0 ) || ( beta > 1.0 ) ||
  806.  
  807.          /* compute alpha & check bounds */
  808.          ( ( alpha = ( u0 - beta * pss->u2 ) * pss->inv_u1 )
  809.             < 0.0 ) ) {
  810.  
  811.         /* whew! missed! */
  812.         goto NextTri ;
  813.         }
  814.     } else {
  815.         /* 0 and 2 vertices have same X value */
  816.  
  817.         /* compute beta and check bounds */
  818.         /* we use "<= 0.0" so that points on the shared interior edge
  819.          * will (generally) be inside only one polygon.
  820.          */
  821.         beta = u0 * pss->inv_u2 ;
  822.         if ( ( beta <= 0.0 ) || ( beta >= 1.0 ) ||
  823.  
  824.          /* compute alpha and check bounds */
  825.          ( ( alpha = ( v0 - beta * pss->v2 ) * pss->inv_v1 )
  826.             < 0.0 ) ) {
  827.  
  828.         /* whew! missed! */
  829.         goto NextTri ;
  830.         }
  831.     }
  832.  
  833.     /* check gamma */
  834.     if ( alpha + beta <= 1.0 ) {
  835.         /* survived */
  836.         inside_flag = !inside_flag ;
  837.     }
  838.  
  839.     NextTri: ;
  840.     }
  841.  
  842.     return( inside_flag ) ;
  843. }
  844.  
  845. void SpackmanCleanup( p_spackman_set )
  846. pSpackmanSet    p_spackman_set ;
  847. {
  848.     free( p_spackman_set ) ;
  849. }
  850.  
  851. /* ======= Trapezoid (bin) algorithm ======================================= */
  852.  
  853. /* Split polygons along set of y bins and sorts the edge fragments.  Testing
  854.  * is done against these fragments.
  855.  *
  856.  * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices, the
  857.  * number of bins desired _bins_, and a pointer to a trapezoid structure
  858.  * _p_trap_set_.
  859.  * Call testing procedure with 2D polygon _pgon_ with _numverts_ number of
  860.  * vertices, _p_trap_set_ pointer to trapezoid structure, and test point
  861.  * _point_, returns 1 if inside, 0 if outside.
  862.  * Call cleanup with pointer to trapezoid structure to free space.
  863.  */
  864. void TrapezoidSetup( pgon, numverts, bins, p_trap_set )
  865. double    pgon[][2] ;
  866. int    numverts ;
  867. int    bins ;
  868. pTrapezoidSet    p_trap_set ;
  869. {
  870. double    *vtx0, *vtx1, *vtxa, *vtxb, slope ;
  871. int    i, j, bin_tot[TOT_BINS], ba, bb, id, full_cross, count ;
  872. double    fba, fbb, vx0, vx1, dy, vy0 ;
  873.  
  874.     p_trap_set->bins = bins ;
  875.     p_trap_set->trapz = (pTrapezoid)malloc( p_trap_set->bins *
  876.         sizeof(Trapezoid));
  877.     MALLOC_CHECK( p_trap_set->trapz ) ;
  878.  
  879.     p_trap_set->minx =
  880.     p_trap_set->maxx = pgon[0][X] ;
  881.     p_trap_set->miny =
  882.     p_trap_set->maxy = pgon[0][Y] ;
  883.  
  884.     for ( i = 1 ; i < numverts ; i++ ) {
  885.     if ( p_trap_set->minx > (vx0 = pgon[i][X]) ) {
  886.         p_trap_set->minx = vx0 ;
  887.     } else if ( p_trap_set->maxx < vx0 ) {
  888.         p_trap_set->maxx = vx0 ;
  889.     }
  890.  
  891.     if ( p_trap_set->miny > (vy0 = pgon[i][Y]) ) {
  892.         p_trap_set->miny = vy0 ;
  893.     } else if ( p_trap_set->maxy < vy0 ) {
  894.         p_trap_set->maxy = vy0 ;
  895.     }
  896.     }
  897.  
  898.     /* add a little to the bounds to ensure everything falls inside area */
  899.     p_trap_set->miny -= EPSILON * (p_trap_set->maxy-p_trap_set->miny) ;
  900.     p_trap_set->maxy += EPSILON * (p_trap_set->maxy-p_trap_set->miny) ;
  901.  
  902.     p_trap_set->ydelta =
  903.         (p_trap_set->maxy-p_trap_set->miny) / (double)p_trap_set->bins ;
  904.     p_trap_set->inv_ydelta = 1.0 / p_trap_set->ydelta ;
  905.  
  906.     /* find how many locations to allocate for each bin */
  907.     for ( i = 0 ; i < p_trap_set->bins ; i++ ) {
  908.     bin_tot[i] = 0 ;
  909.     }
  910.  
  911.     vtx0 = pgon[numverts-1] ;
  912.     for ( i = 0 ; i < numverts ; i++ ) {
  913.     vtx1 = pgon[i] ;
  914.  
  915.     /* skip if Y's identical (edge has no effect) */
  916.     if ( vtx0[Y] != vtx1[Y] ) {
  917.  
  918.         if ( vtx0[Y] < vtx1[Y] ) {
  919.         vtxa = vtx0 ;
  920.         vtxb = vtx1 ;
  921.         } else {
  922.         vtxa = vtx1 ;
  923.         vtxb = vtx0 ;
  924.         }
  925.         ba = (int)(( vtxa[Y]-p_trap_set->miny ) * p_trap_set->inv_ydelta) ;
  926.         fbb = ( vtxb[Y] - p_trap_set->miny ) * p_trap_set->inv_ydelta ;
  927.         bb = (int)fbb ;
  928.         /* if high vertex ends on a boundary, don't go into next boundary */
  929.         if ( fbb - (double)bb == 0.0 ) {
  930.         bb-- ;
  931.         }
  932.  
  933.         /* mark the bins with this edge */
  934.         for ( j = ba ; j <= bb ; j++ ) {
  935.         bin_tot[j]++ ;
  936.         }
  937.     }
  938.  
  939.     vtx0 = vtx1 ;
  940.     }
  941.  
  942.     /* allocate the bin contents and fill in some basics */
  943.     for ( i = 0 ; i < p_trap_set->bins ; i++ ) {
  944.     p_trap_set->trapz[i].edge_set =
  945.         (pEdge*)malloc( bin_tot[i] * sizeof(pEdge) ) ;
  946.     MALLOC_CHECK( p_trap_set->trapz[i].edge_set ) ;
  947.     for ( j = 0 ; j < bin_tot[i] ; j++ ) {
  948.         p_trap_set->trapz[i].edge_set[j] =
  949.         (pEdge)malloc( sizeof(Edge) ) ;
  950.         MALLOC_CHECK( p_trap_set->trapz[i].edge_set[j] ) ;
  951.     }
  952.  
  953.     /* start these off at some awful values; refined below */
  954.     p_trap_set->trapz[i].minx = p_trap_set->maxx ;
  955.     p_trap_set->trapz[i].maxx = p_trap_set->minx ;
  956.     p_trap_set->trapz[i].count = 0 ;
  957.     }
  958.  
  959.     /* now go through list yet again, putting edges in bins */
  960.     vtx0 = pgon[numverts-1] ;
  961.     id = numverts-1 ;
  962.     for ( i = 0 ; i < numverts ; i++ ) {
  963.     vtx1 = pgon[i] ;
  964.  
  965.     /* we can skip edge if Y's are equal */
  966.     if ( vtx0[Y] != vtx1[Y] ) {
  967.         if ( vtx0[Y] < vtx1[Y] ) {
  968.         vtxa = vtx0 ;
  969.         vtxb = vtx1 ;
  970.         } else {
  971.         vtxa = vtx1 ;
  972.         vtxb = vtx0 ;
  973.         }
  974.         fba = ( vtxa[Y] - p_trap_set->miny ) * p_trap_set->inv_ydelta ;
  975.         ba = (int)fba ;
  976.         fbb = ( vtxb[Y] - p_trap_set->miny ) * p_trap_set->inv_ydelta ;
  977.         bb = (int)fbb ;
  978.         /* if high vertex ends on a boundary, don't go into it */
  979.         if ( fbb == (double)bb ) {
  980.         bb-- ;
  981.         }
  982.  
  983.         vx0 = vtxa[X] ;
  984.         dy = vtxa[Y] - vtxb[Y] ;
  985.         slope = p_trap_set->ydelta * ( vtxa[X] - vtxb[X] ) / dy ;
  986.  
  987.         /* set vx1 in case loop is not entered */
  988.         vx1 = vx0 ;
  989.         full_cross = 0 ;
  990.  
  991.         for ( j = ba ; j < bb ; j++, vx0 = vx1 ) {
  992.         /* could increment vx1, but for greater accuracy recompute it */
  993.         vx1 = vtxa[X] + ( (double)(j+1) - fba ) * slope ;
  994.  
  995.         count = p_trap_set->trapz[j].count++ ;
  996.         p_trap_set->trapz[j].edge_set[count]->id = id ;
  997.         p_trap_set->trapz[j].edge_set[count]->full_cross = full_cross;
  998.         TrapBound( j, count, vx0, vx1, p_trap_set ) ;
  999.         full_cross = 1 ;
  1000.         }
  1001.  
  1002.         /* at last bin - fill as above, but with vx1 = vtxb[X] */
  1003.         vx0 = vx1 ;
  1004.         vx1 = vtxb[X] ;
  1005.         count = p_trap_set->trapz[bb].count++ ;
  1006.         p_trap_set->trapz[bb].edge_set[count]->id = id ;
  1007.         /* the last bin is never a full crossing */
  1008.         p_trap_set->trapz[bb].edge_set[count]->full_cross = 0 ;
  1009.         TrapBound( bb, count, vx0, vx1, p_trap_set ) ;
  1010.     }
  1011.  
  1012.     vtx0 = vtx1 ;
  1013.     id = i ;
  1014.     }
  1015.  
  1016.     /* finally, sort the bins' contents by minx */
  1017.     for ( i = 0 ; i < p_trap_set->bins ; i++ ) {
  1018.     qsort( p_trap_set->trapz[i].edge_set, p_trap_set->trapz[i].count,
  1019.         sizeof(pEdge), CompareEdges ) ;
  1020.     }
  1021. }
  1022.  
  1023. void TrapBound( j, count, vx0, vx1, p_trap_set )
  1024. int    j, count ;
  1025. double    vx0, vx1 ;
  1026. pTrapezoidSet    p_trap_set ;
  1027. {
  1028. double    xt ;
  1029.  
  1030.     if ( vx0 > vx1 ) {
  1031.     xt = vx0 ;
  1032.     vx0 = vx1 ;
  1033.     vx1 = xt ;
  1034.     }
  1035.  
  1036.     if ( p_trap_set->trapz[j].minx > vx0 ) {
  1037.     p_trap_set->trapz[j].minx = vx0 ;
  1038.     }
  1039.     if ( p_trap_set->trapz[j].maxx < vx1 ) {
  1040.     p_trap_set->trapz[j].maxx = vx1 ;
  1041.     }
  1042.     p_trap_set->trapz[j].edge_set[count]->minx = vx0 ;
  1043.     p_trap_set->trapz[j].edge_set[count]->maxx = vx1 ;
  1044. }
  1045.  
  1046. /* used by qsort to sort */
  1047. int CompareEdges( u, v )
  1048. pEdge    *u, *v ;
  1049. {
  1050.     if ( (*u)->minx == (*v)->minx ) {
  1051.     return( 0 ) ;
  1052.     } else {
  1053.     return( (*u)->minx < (*v)->minx ? -1 : 1 ) ;
  1054.     }
  1055. }
  1056.  
  1057. int TrapezoidTest( pgon, numverts, p_trap_set, point )
  1058. double    pgon[][2] ;
  1059. int    numverts ;
  1060. pTrapezoidSet    p_trap_set ;
  1061. double    point[2] ;
  1062. {
  1063. int    j, b, count, id ;
  1064. double    tx, ty, *vtx0, *vtx1 ;
  1065. pEdge    *pp_bin ;
  1066. pTrapezoid    p_trap ;
  1067. int    inside_flag ;
  1068.  
  1069.     inside_flag = 0 ;
  1070.  
  1071.     /* first, is point inside bounding rectangle? */
  1072.     if ( ( ty = point[Y] ) < p_trap_set->miny ||
  1073.      ty >= p_trap_set->maxy ||
  1074.      ( tx = point[X] ) < p_trap_set->minx ||
  1075.      tx >= p_trap_set->maxx ) {
  1076.  
  1077.     /* outside of box */
  1078.     return( 0 ) ;
  1079.     }
  1080.  
  1081.     /* what bin are we in? */
  1082.     b = ( ty - p_trap_set->miny ) * p_trap_set->inv_ydelta ;
  1083.  
  1084.     /* find if we're inside this bin's bounds */
  1085.     if ( tx < (p_trap = &p_trap_set->trapz[b])->minx ||
  1086.      tx > p_trap->maxx ) {
  1087.  
  1088.     /* outside of box */
  1089.     return( 0 ) ;
  1090.     }
  1091.  
  1092.     /* now search bin for crossings */
  1093.     pp_bin = p_trap->edge_set ;
  1094.     count = p_trap->count ;
  1095.     for ( j = 0 ; j < count ; j++, pp_bin++ ) {
  1096.     if ( tx < (*pp_bin)->minx ) {
  1097.  
  1098.         /* all remaining edges are to right of point, so test them */
  1099.         do {
  1100.         if ( (*pp_bin)->full_cross ) {
  1101.             inside_flag = !inside_flag ;
  1102.         } else {
  1103.             id = (*pp_bin)->id ;
  1104.             if ( ( ty <= pgon[id][Y] ) !=
  1105.                 ( ty <= pgon[(id+1)%numverts][Y] ) ) {
  1106.  
  1107.             /* point crosses edge in Y, so must cross */
  1108.             inside_flag = !inside_flag ;
  1109.             }
  1110.         }
  1111.         pp_bin++ ;
  1112.         } while ( ++j < count ) ;
  1113.         goto Exit;
  1114.  
  1115.     } else if ( tx < (*pp_bin)->maxx ) {
  1116.         /* edge is overlapping point in X, check it */
  1117.         id = (*pp_bin)->id ;
  1118.         vtx0 = pgon[id] ;
  1119.         vtx1 = pgon[(id+1)%numverts] ;
  1120.  
  1121.         if ( (*pp_bin)->full_cross ||
  1122.          ( ty <= vtx0[Y] ) != ( ty <= vtx1[Y] ) ) {
  1123.  
  1124.         /* edge crosses in Y, so have to do full crossings test */
  1125.         if ( (vtx0[X] -
  1126.             (vtx0[Y] - ty )*
  1127.             ( vtx1[X]-vtx0[X])/(vtx1[Y]-vtx0[Y])) >= tx ) {
  1128.             inside_flag = !inside_flag ;
  1129.         }
  1130.         }
  1131.  
  1132.     } /* else edge is to left of point, ignore it */
  1133.     }
  1134.  
  1135.     Exit:
  1136.     return( inside_flag ) ;
  1137. }
  1138.  
  1139. void TrapezoidCleanup( p_trap_set )
  1140. pTrapezoidSet    p_trap_set ;
  1141. {
  1142. int    i, j, count ;
  1143.  
  1144.     for ( i = 0 ; i < p_trap_set->bins ; i++ ) {
  1145.     /* all of these should have bin sets, but check just in case */
  1146.     if ( p_trap_set->trapz[i].edge_set ) {
  1147.         count = p_trap_set->trapz[i].count ;
  1148.         for ( j = 0 ; j < count ; j++ ) {
  1149.         if ( p_trap_set->trapz[i].edge_set[j] ) {
  1150.             free( p_trap_set->trapz[i].edge_set[j] ) ;
  1151.         }
  1152.         }
  1153.         free( p_trap_set->trapz[i].edge_set ) ;
  1154.     }
  1155.     }
  1156.     free( p_trap_set->trapz ) ;
  1157. }
  1158.  
  1159. /* ======= Grid algorithm ================================================= */
  1160.  
  1161. /* Impose a grid upon the polygon and test only the local edges against the
  1162.  * point.
  1163.  *
  1164.  * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
  1165.  * grid resolution _resolution_ and a pointer to a grid structure _p_gs_.
  1166.  * Call testing procedure with a pointer to this array and test point _point_,
  1167.  * returns 1 if inside, 0 if outside.
  1168.  * Call cleanup with pointer to grid structure to free space.
  1169.  */
  1170.  
  1171. /* Strategy for setup:
  1172.  *   Get bounds of polygon, allocate grid.
  1173.  *   "Walk" each edge of the polygon and note which edges have been crossed
  1174.  *     and what cells are entered (points on a grid edge are always considered
  1175.  *     to be above that edge).    Keep a record of the edges overlapping a cell.
  1176.  *     For cells with edges, determine if any cell border has no edges passing
  1177.  *     through it and so can be used for shooting a test ray.
  1178.  *     Keep track of the parity of the x (horizontal) grid cell borders for
  1179.  *     use in determining whether the grid corners are inside or outside.
  1180.  */
  1181. void GridSetup( pgon, numverts, resolution, p_gs )
  1182. double    pgon[][2] ;
  1183. int    numverts ;
  1184. int    resolution ;
  1185. pGridSet    p_gs ;
  1186. {
  1187. double    *vtx0, *vtx1, *vtxa, *vtxb, *p_gl ;
  1188. int    i, j, gc_clear_flags ;
  1189. double    vx0, vx1, vy0, vy1, gxdiff, gydiff, eps ;
  1190. pGridCell    p_gc, p_ngc ;
  1191. double    xdiff, ydiff, tmax, inv_x, inv_y, xdir, ydir, t_near, tx, ty ;
  1192. double    tgcx, tgcy ;
  1193. int    gcx, gcy, sign_x ;
  1194. int    y_flag, io_state ;
  1195.  
  1196.     p_gs->xres = p_gs->yres = resolution ;
  1197.     p_gs->tot_cells = p_gs->xres * p_gs->yres ;
  1198.     p_gs->glx = (double *)malloc( (p_gs->xres+1) * sizeof(double));
  1199.     MALLOC_CHECK( p_gs->glx ) ;
  1200.     p_gs->gly = (double *)malloc( (p_gs->yres+1) * sizeof(double));
  1201.     MALLOC_CHECK( p_gs->gly ) ;
  1202.     p_gs->gc = (pGridCell)malloc( p_gs->tot_cells * sizeof(GridCell));
  1203.     MALLOC_CHECK( p_gs->gc ) ;
  1204.  
  1205.     p_gs->minx =
  1206.     p_gs->maxx = pgon[0][X] ;
  1207.     p_gs->miny =
  1208.     p_gs->maxy = pgon[0][Y] ;
  1209.  
  1210.     /* find bounds of polygon */
  1211.     for ( i = 1 ; i < numverts ; i++ ) {
  1212.     vx0 = pgon[i][X] ;
  1213.     if ( p_gs->minx > vx0 ) {
  1214.         p_gs->minx = vx0 ;
  1215.     } else if ( p_gs->maxx < vx0 ) {
  1216.         p_gs->maxx = vx0 ;
  1217.     }
  1218.  
  1219.     vy0 = pgon[i][Y] ;
  1220.     if ( p_gs->miny > vy0 ) {
  1221.         p_gs->miny = vy0 ;
  1222.     } else if ( p_gs->maxy < vy0 ) {
  1223.         p_gs->maxy = vy0 ;
  1224.     }
  1225.     }
  1226.  
  1227.     /* add a little to the bounds to ensure everything falls inside area */
  1228.     gxdiff = p_gs->maxx - p_gs->minx ;
  1229.     gydiff = p_gs->maxy - p_gs->miny ;
  1230.     p_gs->minx -= EPSILON * gxdiff ;
  1231.     p_gs->maxx += EPSILON * gxdiff ;
  1232.     p_gs->miny -= EPSILON * gydiff ;
  1233.     p_gs->maxy += EPSILON * gydiff ;
  1234.  
  1235.     /* avoid roundoff problems near corners by not getting too close to them */
  1236.     eps = 1e-9 * ( gxdiff + gydiff ) ;
  1237.  
  1238.     /* use the new bounds to compute cell widths */
  1239.     TryAgain:
  1240.     p_gs->xdelta =
  1241.         (p_gs->maxx-p_gs->minx) / (double)p_gs->xres ;
  1242.     p_gs->inv_xdelta = 1.0 / p_gs->xdelta ;
  1243.  
  1244.     p_gs->ydelta =
  1245.         (p_gs->maxy-p_gs->miny) / (double)p_gs->yres ;
  1246.     p_gs->inv_ydelta = 1.0 / p_gs->ydelta ;
  1247.  
  1248.     for ( i = 0, p_gl = p_gs->glx ; i < p_gs->xres ; i++ ) {
  1249.     *p_gl++ = p_gs->minx + i * p_gs->xdelta ;
  1250.     }
  1251.     /* make last grid corner precisely correct */
  1252.     *p_gl = p_gs->maxx ;
  1253.  
  1254.     for ( i = 0, p_gl = p_gs->gly ; i < p_gs->yres ; i++ ) {
  1255.     *p_gl++ = p_gs->miny + i * p_gs->ydelta ;
  1256.     }
  1257.     *p_gl = p_gs->maxy ;
  1258.  
  1259.     for ( i = 0, p_gc = p_gs->gc ; i < p_gs->tot_cells ; i++, p_gc++ ) {
  1260.     p_gc->tot_edges = 0 ;
  1261.     p_gc->gc_flags = 0x0 ;
  1262.     p_gc->gr = NULL ;
  1263.     }
  1264.  
  1265.     /* loop through edges and insert into grid structure */
  1266.     vtx0 = pgon[numverts-1] ;
  1267.     for ( i = 0 ; i < numverts ; i++ ) {
  1268.     vtx1 = pgon[i] ;
  1269.  
  1270.     if ( vtx0[Y] < vtx1[Y] ) {
  1271.         vtxa = vtx0 ;
  1272.         vtxb = vtx1 ;
  1273.     } else {
  1274.         vtxa = vtx1 ;
  1275.         vtxb = vtx0 ;
  1276.     }
  1277.  
  1278.     /* Set x variable for the direction of the ray */
  1279.     xdiff = vtxb[X] - vtxa[X] ;
  1280.     ydiff = vtxb[Y] - vtxa[Y] ;
  1281.     tmax = sqrt( xdiff * xdiff + ydiff * ydiff ) ;
  1282.  
  1283.     /* if edge is of 0 length, ignore it (useless edge) */
  1284.     if ( tmax == 0.0 ) goto NextEdge ;
  1285.  
  1286.     inv_y = tmax / ydiff ;
  1287.     xdir = xdiff / tmax ;
  1288.     ydir = ydiff / tmax ;
  1289.  
  1290.     gcx = (int)(( vtxa[X] - p_gs->minx ) * p_gs->inv_xdelta) ;
  1291.     gcy = (int)(( vtxa[Y] - p_gs->miny ) * p_gs->inv_ydelta) ;
  1292.  
  1293.     /* get information about slopes of edge, etc */
  1294.     if ( vtxa[X] == vtxb[X] ) {
  1295.         sign_x = 0 ;
  1296.         tx = HUGE ;
  1297.     } else {
  1298.         inv_x = tmax / xdiff ;
  1299.         tx = p_gs->xdelta * (double)gcx + p_gs->minx - vtxa[X] ;
  1300.         if ( vtxa[X] < vtxb[X] ) {
  1301.         sign_x = 1 ;
  1302.         tx += p_gs->xdelta ;
  1303.         tgcx = p_gs->xdelta * inv_x ;
  1304.         } else {
  1305.         sign_x = -1 ;
  1306.         tgcx = -p_gs->xdelta * inv_x ;
  1307.         }
  1308.         tx *= inv_x ;
  1309.     }
  1310.  
  1311.     if ( vtxa[Y] == vtxb[Y] ) {
  1312.         ty = HUGE ;
  1313.     } else {
  1314.         inv_y = tmax / ydiff ;
  1315.         ty = (p_gs->ydelta * (double)(gcy+1) + p_gs->miny - vtxa[Y])
  1316.         * inv_y ;
  1317.         tgcy = p_gs->ydelta * inv_y ;
  1318.     }
  1319.  
  1320.     p_gc = &p_gs->gc[gcy*p_gs->xres+gcx] ;
  1321.  
  1322.     vx0 = vtxa[X] ;
  1323.     vy0 = vtxa[Y] ;
  1324.  
  1325.     t_near = 0.0 ;
  1326.  
  1327.     do {
  1328.         /* choose the next boundary, but don't move yet */
  1329.         if ( tx <= ty ) {
  1330.         gcx += sign_x ;
  1331.  
  1332.         ty -= tx ;
  1333.         t_near += tx ;
  1334.         tx = tgcx ;
  1335.  
  1336.         /* note which edge is hit when leaving this cell */
  1337.         if ( t_near < tmax ) {
  1338.             if ( sign_x > 0 ) {
  1339.             p_gc->gc_flags |= GC_R_EDGE_HIT ;
  1340.             vx1 = p_gs->glx[gcx] ;
  1341.             } else {
  1342.             p_gc->gc_flags |= GC_L_EDGE_HIT ;
  1343.             vx1 = p_gs->glx[gcx+1] ;
  1344.             }
  1345.  
  1346.             /* get new location */
  1347.             vy1 = t_near * ydir + vtxa[Y] ;
  1348.         } else {
  1349.             /* end of edge, so get exact value */
  1350.             vx1 = vtxb[X] ;
  1351.             vy1 = vtxb[Y] ;
  1352.         }
  1353.  
  1354.         y_flag = FALSE ;
  1355.  
  1356.         } else {
  1357.  
  1358.         gcy++ ;
  1359.  
  1360.         tx -= ty ;
  1361.         t_near += ty ;
  1362.         ty = tgcy ;
  1363.  
  1364.         /* note top edge is hit when leaving this cell */
  1365.         if ( t_near < tmax ) {
  1366.             p_gc->gc_flags |= GC_T_EDGE_HIT ;
  1367.             /* this toggles the parity bit */
  1368.             p_gc->gc_flags ^= GC_T_EDGE_PARITY ;
  1369.  
  1370.             /* get new location */
  1371.             vx1 = t_near * xdir + vtxa[X] ;
  1372.             vy1 = p_gs->gly[gcy] ;
  1373.         } else {
  1374.             /* end of edge, so get exact value */
  1375.             vx1 = vtxb[X] ;
  1376.             vy1 = vtxb[Y] ;
  1377.         }
  1378.  
  1379.         y_flag = TRUE ;
  1380.         }
  1381.  
  1382.         /* check for corner crossing, then mark the cell we're in */
  1383.         if ( !AddGridRecAlloc( p_gc, vx0, vy0, vx1, vy1, eps ) ) {
  1384.         /* warning, danger - we have just crossed a corner.
  1385.          * There are all kinds of topological messiness we could
  1386.          * do to get around this case, but they're a headache.
  1387.          * The simplest recovery is just to change the extents a bit
  1388.          * and redo the meshing, so that hopefully no edges will
  1389.          * perfectly cross a corner.  Since it's a preprocess, we
  1390.          * don't care too much about the time to do it.
  1391.          */
  1392.  
  1393.         /* clean out all grid records */
  1394.         for ( i = 0, p_gc = p_gs->gc
  1395.             ; i < p_gs->tot_cells
  1396.             ; i++, p_gc++ ) {
  1397.  
  1398.             if ( p_gc->gr ) {
  1399.             free( p_gc->gr ) ;
  1400.             }
  1401.         }
  1402.  
  1403.         /* make the bounding box ever so slightly larger, hopefully
  1404.          * changing the alignment of the corners.
  1405.          */
  1406.         p_gs->minx -= EPSILON * gxdiff * 0.24 ;
  1407.         p_gs->miny -= EPSILON * gydiff * 0.10 ;
  1408.  
  1409.         /* yes, it's the dreaded goto - run in fear for your lives! */
  1410.         goto TryAgain ;
  1411.         }
  1412.  
  1413.         if ( t_near < tmax ) {
  1414.         /* note how we're entering the next cell */
  1415.         /* TBD: could be done faster by incrementing index in the
  1416.          * incrementing code, above */
  1417.         p_gc = &p_gs->gc[gcy*p_gs->xres+gcx] ;
  1418.  
  1419.         if ( y_flag ) {
  1420.             p_gc->gc_flags |= GC_B_EDGE_HIT ;
  1421.             /* this toggles the parity bit */
  1422.             p_gc->gc_flags ^= GC_B_EDGE_PARITY ;
  1423.         } else {
  1424.             p_gc->gc_flags |=
  1425.             ( sign_x > 0 ) ? GC_L_EDGE_HIT : GC_R_EDGE_HIT ;
  1426.         }
  1427.         }
  1428.  
  1429.         vx0 = vx1 ;
  1430.         vy0 = vy1 ;
  1431.     }
  1432.     /* have we gone further than the end of the edge? */
  1433.     while ( t_near < tmax ) ;
  1434.  
  1435.     NextEdge:
  1436.     vtx0 = vtx1 ;
  1437.     }
  1438.  
  1439.     /* the grid is all set up, now set up the inside/outside value of each
  1440.      * corner.
  1441.      */
  1442.     p_gc = p_gs->gc ;
  1443.     p_ngc = &p_gs->gc[p_gs->xres] ;
  1444.  
  1445.     /* we know the bottom and top rows are all outside, so no flag is set */
  1446.     for ( i = 1; i < p_gs->yres ; i++ ) {
  1447.     /* start outside */
  1448.     io_state = 0x0 ;
  1449.  
  1450.     for ( j = 0; j < p_gs->xres ; j++ ) {
  1451.  
  1452.         if ( io_state ) {
  1453.         /* change cell left corners to inside */
  1454.         p_gc->gc_flags |= GC_TL_IN ;
  1455.         p_ngc->gc_flags |= GC_BL_IN ;
  1456.         }
  1457.  
  1458.         if ( p_gc->gc_flags & GC_T_EDGE_PARITY ) {
  1459.         io_state = !io_state ;
  1460.         }
  1461.  
  1462.         if ( io_state ) {
  1463.         /* change cell right corners to inside */
  1464.         p_gc->gc_flags |= GC_TR_IN ;
  1465.         p_ngc->gc_flags |= GC_BR_IN ;
  1466.         }
  1467.  
  1468.         p_gc++ ;
  1469.         p_ngc++ ;
  1470.     }
  1471.     }
  1472.  
  1473.     p_gc = p_gs->gc ;
  1474.     for ( i = 0; i < p_gs->tot_cells ; i++ ) {
  1475.  
  1476.     /* reverse parity of edge clear (1==edge clear) */
  1477.     gc_clear_flags = p_gc->gc_flags ^ GC_ALL_EDGE_CLEAR ;
  1478.     if ( gc_clear_flags & GC_L_EDGE_CLEAR ) {
  1479.         p_gc->gc_flags |= GC_AIM_L ;
  1480.     } else
  1481.     if ( gc_clear_flags & GC_B_EDGE_CLEAR ) {
  1482.         p_gc->gc_flags |= GC_AIM_B ;
  1483.     } else
  1484.     if ( gc_clear_flags & GC_R_EDGE_CLEAR ) {
  1485.         p_gc->gc_flags |= GC_AIM_R ;
  1486.     } else
  1487.     if ( gc_clear_flags & GC_T_EDGE_CLEAR ) {
  1488.         p_gc->gc_flags |= GC_AIM_T ;
  1489.     } else {
  1490.         /* all edges have something on them, do full test */
  1491.         p_gc->gc_flags |= GC_AIM_C ;
  1492.     }
  1493.     p_gc++ ;
  1494.     }
  1495. }
  1496.  
  1497. int AddGridRecAlloc( p_gc, xa, ya, xb, yb, eps )
  1498. pGridCell    p_gc ;
  1499. double    xa,ya,xb,yb,eps ;
  1500. {
  1501. pGridRec    p_gr ;
  1502. double        slope, inv_slope ;
  1503.  
  1504.     if ( Near(ya, yb, eps) ) {
  1505.     if ( Near(xa, xb, eps) ) {
  1506.         /* edge is 0 length, so get rid of it */
  1507.         return( FALSE ) ;
  1508.     } else {
  1509.         /* horizontal line */
  1510.         slope = HUGE ;
  1511.         inv_slope = 0.0 ;
  1512.     }
  1513.     } else {
  1514.     if ( Near(xa, xb, eps) ) {
  1515.         /* vertical line */
  1516.         slope = 0.0 ;
  1517.         inv_slope = HUGE ;
  1518.     } else {
  1519.         slope = (xb-xa)/(yb-ya) ;
  1520.         inv_slope = (yb-ya)/(xb-xa) ;
  1521.     }
  1522.     }
  1523.  
  1524.     p_gc->tot_edges++ ;
  1525.     if ( p_gc->tot_edges <= 1 ) {
  1526.     p_gc->gr = (pGridRec)malloc( sizeof(GridRec) ) ;
  1527.     } else {
  1528.     p_gc->gr = (pGridRec)realloc( p_gc->gr,
  1529.         p_gc->tot_edges * sizeof(GridRec) ) ;
  1530.     }
  1531.     MALLOC_CHECK( p_gc->gr ) ;
  1532.     p_gr = &p_gc->gr[p_gc->tot_edges-1] ;
  1533.  
  1534.     p_gr->slope = slope ;
  1535.     p_gr->inv_slope = inv_slope ;
  1536.  
  1537.     p_gr->xa = xa ;
  1538.     p_gr->ya = ya ;
  1539.     if ( xa <= xb ) {
  1540.     p_gr->minx = xa ;
  1541.     p_gr->maxx = xb ;
  1542.     } else {
  1543.     p_gr->minx = xb ;
  1544.     p_gr->maxx = xa ;
  1545.     }
  1546.     if ( ya <= yb ) {
  1547.     p_gr->miny = ya ;
  1548.     p_gr->maxy = yb ;
  1549.     } else {
  1550.     p_gr->miny = yb ;
  1551.     p_gr->maxy = ya ;
  1552.     }
  1553.  
  1554.     /* P2 - P1 */
  1555.     p_gr->ax = xb - xa ;
  1556.     p_gr->ay = yb - ya ;
  1557.  
  1558.     return( TRUE ) ;
  1559. }
  1560.  
  1561. /* Test point against grid and edges in the cell (if any).  Algorithm:
  1562.  *    Check bounding box; if outside then return.
  1563.  *    Check cell point is inside; if simple inside or outside then return.
  1564.  *    Find which edge or corner is considered to be the best for testing and
  1565.  *      send a test ray towards it, counting the crossings.  Add in the
  1566.  *      state of the edge or corner the ray went to and so determine the
  1567.  *      state of the point (inside or outside).
  1568.  */
  1569. int GridTest( p_gs, point )
  1570. register pGridSet    p_gs ;
  1571. double    point[2] ;
  1572. {
  1573. int    j, count, init_flag ;
  1574. pGridCell    p_gc ;
  1575. pGridRec    p_gr ;
  1576. double    tx, ty, xcell, ycell, bx,by,cx,cy, cornerx, cornery ;
  1577. double    alpha, beta, denom ;
  1578. unsigned short    gc_flags ;
  1579. int    inside_flag ;
  1580.  
  1581.     /* first, is point inside bounding rectangle? */
  1582.     if ( ( ty = point[Y] ) < p_gs->miny ||
  1583.      ty >= p_gs->maxy ||
  1584.      ( tx = point[X] ) < p_gs->minx ||
  1585.      tx >= p_gs->maxx ) {
  1586.  
  1587.     /* outside of box */
  1588.     inside_flag = FALSE ;
  1589.     } else {
  1590.  
  1591.     /* what cell are we in? */
  1592.     ycell = ( ty - p_gs->miny ) * p_gs->inv_ydelta ;
  1593.     xcell = ( tx - p_gs->minx ) * p_gs->inv_xdelta ;
  1594.     p_gc = &p_gs->gc[((int)ycell)*p_gs->xres + (int)xcell] ;
  1595.  
  1596.     /* is cell simple? */
  1597.     count = p_gc->tot_edges ;
  1598.     if ( count ) {
  1599.         /* no, so find an edge which is free. */
  1600.         gc_flags = p_gc->gc_flags ;
  1601.         p_gr = p_gc->gr ;
  1602.         switch( gc_flags & GC_AIM ) {
  1603.         case GC_AIM_L:
  1604.         /* left edge is clear, shoot X- ray */
  1605.         /* note - this next statement requires that GC_BL_IN is 1 */
  1606.         inside_flag = gc_flags & GC_BL_IN ;
  1607.         for ( j = count+1 ; --j ; p_gr++ ) {
  1608.             /* test if y is between edges */
  1609.             if ( ty >= p_gr->miny && ty < p_gr->maxy ) {
  1610.             if ( tx > p_gr->maxx ) {
  1611.                 inside_flag = !inside_flag ;
  1612.             } else if ( tx > p_gr->minx ) {
  1613.                 /* full computation */
  1614.                 if ( ( p_gr->xa -
  1615.                 ( p_gr->ya - ty ) * p_gr->slope ) < tx ) {
  1616.                 inside_flag = !inside_flag ;
  1617.                 }
  1618.             }
  1619.             }
  1620.         }
  1621.         break ;
  1622.  
  1623.         case GC_AIM_B:
  1624.         /* bottom edge is clear, shoot Y+ ray */
  1625.         /* note - this next statement requires that GC_BL_IN is 1 */
  1626.         inside_flag = gc_flags & GC_BL_IN ;
  1627.         for ( j = count+1 ; --j ; p_gr++ ) {
  1628.             /* test if x is between edges */
  1629.             if ( tx >= p_gr->minx && tx < p_gr->maxx ) {
  1630.             if ( ty > p_gr->maxy ) {
  1631.                 inside_flag = !inside_flag ;
  1632.             } else if ( ty > p_gr->miny ) {
  1633.                 /* full computation */
  1634.                 if ( ( p_gr->ya - ( p_gr->xa - tx ) *
  1635.                     p_gr->inv_slope ) < ty ) {
  1636.                 inside_flag = !inside_flag ;
  1637.                 }
  1638.             }
  1639.             }
  1640.         }
  1641.         break ;
  1642.  
  1643.         case GC_AIM_R:
  1644.         /* right edge is clear, shoot X+ ray */
  1645.         inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0 ;
  1646.  
  1647.         /* TBD: Note, we could have sorted the edges to be tested
  1648.          * by miny or somesuch, and so be able to cut testing
  1649.          * short when the list's miny > point.y .
  1650.          */
  1651.         for ( j = count+1 ; --j ; p_gr++ ) {
  1652.             /* test if y is between edges */
  1653.             if ( ty >= p_gr->miny && ty < p_gr->maxy ) {
  1654.             if ( tx <= p_gr->minx ) {
  1655.                 inside_flag = !inside_flag ;
  1656.             } else if ( tx <= p_gr->maxx ) {
  1657.                 /* full computation */
  1658.                 if ( ( p_gr->xa -
  1659.                 ( p_gr->ya - ty ) * p_gr->slope ) >= tx ) {
  1660.                 inside_flag = !inside_flag ;
  1661.                 }
  1662.             }
  1663.             }
  1664.         }
  1665.         break ;
  1666.  
  1667.         case GC_AIM_T:
  1668.         /* top edge is clear, shoot Y+ ray */
  1669.         inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0 ;
  1670.         for ( j = count+1 ; --j ; p_gr++ ) {
  1671.             /* test if x is between edges */
  1672.             if ( tx >= p_gr->minx && tx < p_gr->maxx ) {
  1673.             if ( ty <= p_gr->miny ) {
  1674.                 inside_flag = !inside_flag ;
  1675.             } else if ( ty <= p_gr->maxy ) {
  1676.                 /* full computation */
  1677.                 if ( ( p_gr->ya - ( p_gr->xa - tx ) *
  1678.                     p_gr->inv_slope ) >= ty ) {
  1679.                 inside_flag = !inside_flag ;
  1680.                 }
  1681.             }
  1682.             }
  1683.         }
  1684.         break ;
  1685.  
  1686.         case GC_AIM_C:
  1687.         /* no edge is clear, bite the bullet and test
  1688.          * against the bottom left corner.
  1689.          * We use Franklin Antonio's algorithm (Graphics Gems III).
  1690.          */
  1691.         /* TBD: Faster yet might be to test against the closest
  1692.          * corner to the cell location, but our hope is that we
  1693.          * rarely need to do this testing at all.
  1694.          */
  1695.         inside_flag = ((gc_flags & GC_BL_IN) == GC_BL_IN) ;
  1696.         init_flag = TRUE ;
  1697.  
  1698.         /* get lower left corner coordinate */
  1699.         cornerx = p_gs->glx[(int)xcell] ;
  1700.         cornery = p_gs->gly[(int)ycell] ;
  1701.         for ( j = count+1 ; --j ; p_gr++ ) {
  1702.  
  1703.             /* quick out test: if test point is
  1704.              * less than minx & miny, edge cannot overlap.
  1705.              */
  1706.             if ( tx >= p_gr->minx && ty >= p_gr->miny ) {
  1707.  
  1708.             /* quick test failed, now check if test point and
  1709.              * corner are on different sides of edge.
  1710.              */
  1711.             if ( init_flag ) {
  1712.                 /* Compute these at most once for test */
  1713.                 /* P3 - P4 */
  1714.                 bx = tx - cornerx ;
  1715.                 by = ty - cornery ;
  1716.                 init_flag = FALSE ;
  1717.             }
  1718.             denom = p_gr->ay * bx - p_gr->ax * by ;
  1719.             if ( denom != 0.0 ) {
  1720.                 /* lines are not collinear, so continue */
  1721.                 /* P1 - P3 */
  1722.                 cx = p_gr->xa - tx ;
  1723.                 cy = p_gr->ya - ty ;
  1724.                 alpha = by * cx - bx * cy ;
  1725.                 if ( denom > 0.0 ) {
  1726.                 if ( alpha < 0.0 || alpha >= denom ) {
  1727.                     /* test edge not hit */
  1728.                     goto NextEdge ;
  1729.                 }
  1730.                 beta = p_gr->ax * cy - p_gr->ay * cx ;
  1731.                 if ( beta < 0.0 || beta >= denom ) {
  1732.                     /* polygon edge not hit */
  1733.                     goto NextEdge ;
  1734.                 }
  1735.                 } else {
  1736.                 if ( alpha > 0.0 || alpha <= denom ) {
  1737.                     /* test edge not hit */
  1738.                     goto NextEdge ;
  1739.                 }
  1740.                 beta = p_gr->ax * cy - p_gr->ay * cx ;
  1741.                 if ( beta > 0.0 || beta <= denom ) {
  1742.                     /* polygon edge not hit */
  1743.                     goto NextEdge ;
  1744.                 }
  1745.                 }
  1746.                 inside_flag = !inside_flag ;
  1747.             }
  1748.  
  1749.             }
  1750.             NextEdge: ;
  1751.         }
  1752.         break ;
  1753.         }
  1754.     } else {
  1755.         /* simple cell, so if lower left corner is in,
  1756.          * then cell is inside.
  1757.          */
  1758.         inside_flag = p_gc->gc_flags & GC_BL_IN ;
  1759.     }
  1760.     }
  1761.  
  1762.     return( inside_flag ) ;
  1763. }
  1764.  
  1765. void GridCleanup( p_gs )
  1766. pGridSet    p_gs ;
  1767. {
  1768. int    i ;
  1769. pGridCell    p_gc ;
  1770.  
  1771.     for ( i = 0, p_gc = p_gs->gc
  1772.     ; i < p_gs->tot_cells
  1773.     ; i++, p_gc++ ) {
  1774.  
  1775.     if ( p_gc->gr ) {
  1776.         free( p_gc->gr ) ;
  1777.     }
  1778.     }
  1779.     free( p_gs->glx ) ;
  1780.     free( p_gs->gly ) ;
  1781.     free( p_gs->gc ) ;
  1782. }
  1783.  
  1784. /* ======= Exterior (convex only) algorithm =============================== */
  1785.  
  1786. /* Test the edges of the convex polygon against the point.  If the point is
  1787.  * outside any edge, the point is outside the polygon.
  1788.  *
  1789.  * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
  1790.  * which returns a pointer to a plane set array.
  1791.  * Call testing procedure with a pointer to this array, _numverts_, and
  1792.  * test point _point_, returns 1 if inside, 0 if outside.
  1793.  * Call cleanup with pointer to plane set array to free space.
  1794.  *
  1795.  * RANDOM can be defined for this test.
  1796.  * CONVEX must be defined for this test; it is not usable for general polygons.
  1797.  */
  1798.  
  1799. #ifdef    CONVEX
  1800. /* make exterior plane set */
  1801. pPlaneSet ExteriorSetup( pgon, numverts )
  1802. double    pgon[][2] ;
  1803. int    numverts ;
  1804. {
  1805. int    p1, p2, flip_edge ;
  1806. pPlaneSet    pps, pps_return ;
  1807. #ifdef    RANDOM
  1808. int    i, ind ;
  1809. PlaneSet    ps_temp ;
  1810. #endif
  1811.  
  1812.     pps = pps_return =
  1813.         (pPlaneSet)malloc( numverts * sizeof( PlaneSet )) ;
  1814.     MALLOC_CHECK( pps ) ;
  1815.  
  1816.     /* take cross product of vertex to find handedness */
  1817.     flip_edge = (pgon[0][X] - pgon[1][X]) * (pgon[1][Y] - pgon[2][Y] ) >
  1818.         (pgon[0][Y] - pgon[1][Y]) * (pgon[1][X] - pgon[2][X] ) ;
  1819.  
  1820.     /* Generate half-plane boundary equations now for faster testing later.
  1821.      * vx & vy are the edge's normal, c is the offset from the origin.
  1822.      */
  1823.     for ( p1 = numverts-1, p2 = 0 ; p2 < numverts ; p1 = p2, p2++, pps++ ) {
  1824.     pps->vx = pgon[p1][Y] - pgon[p2][Y] ;
  1825.     pps->vy = pgon[p2][X] - pgon[p1][X] ;
  1826.     pps->c = pps->vx * pgon[p1][X] + pps->vy * pgon[p1][Y] ;
  1827.  
  1828.     /* check sense and reverse plane edge if need be */
  1829.     if ( flip_edge ) {
  1830.         pps->vx = -pps->vx ;
  1831.         pps->vy = -pps->vy ;
  1832.         pps->c  = -pps->c ;
  1833.     }
  1834.     }
  1835.  
  1836. #ifdef    RANDOM
  1837.     /* Randomize the order of the edges to improve chance of early out */
  1838.     /* There are better orders, but the default order is the worst */
  1839.     for ( i = 0, pps = pps_return
  1840.     ; i < numverts
  1841.     ; i++ ) {
  1842.  
  1843.     ind = (int)(RAN01() * numverts ) ;
  1844.     if ( ( ind < 0 ) || ( ind >= numverts ) ) {
  1845.         fprintf( stderr,
  1846.         "Yikes, the random number generator is returning values\n" ) ;
  1847.         fprintf( stderr,
  1848.         "outside the range [0.0,1.0), so please fix the code!\n" ) ;
  1849.         ind = 0 ;
  1850.     }
  1851.  
  1852.     /* swap edges */
  1853.     ps_temp = *pps ;
  1854.     *pps = pps_return[ind] ;
  1855.     pps_return[ind] = ps_temp ;
  1856.     }
  1857. #endif
  1858.  
  1859.     return( pps_return ) ;
  1860. }
  1861.  
  1862. /* Check point for outside of all planes */
  1863. /* note that we don't need "pgon", since it's been processed into
  1864.  * its corresponding PlaneSet.
  1865.  */
  1866. int ExteriorTest( p_ext_set, numverts, point )
  1867. pPlaneSet    p_ext_set ;
  1868. int    numverts ;
  1869. double    point[2] ;
  1870. {
  1871. register PlaneSet    *pps ;
  1872. register int    p0 ;
  1873. register double tx, ty ;
  1874. int    inside_flag ;
  1875.  
  1876.     tx = point[X] ;
  1877.     ty = point[Y] ;
  1878.  
  1879.     for ( p0 = numverts+1, pps = p_ext_set ; --p0 ; pps++ ) {
  1880.  
  1881.     /* test if the point is outside this edge */
  1882.     if ( pps->vx * tx + pps->vy * ty > pps->c ) {
  1883.         return( 0 ) ;
  1884.     }
  1885.     }
  1886.     /* if we make it to here, we were inside all edges */
  1887.     return( 1 ) ;
  1888. }
  1889.  
  1890. void ExteriorCleanup( p_ext_set )
  1891. pPlaneSet    p_ext_set ;
  1892. {
  1893.     free( p_ext_set ) ;
  1894. }
  1895. #endif
  1896.  
  1897. /* ======= Inclusion (convex only) algorithm ============================== */
  1898.  
  1899. /* Create an efficiency structure (see Preparata) for the convex polygon which
  1900.  * allows binary searching to find which edge to test the point against.  This
  1901.  * algorithm is O(log n).
  1902.  *
  1903.  * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
  1904.  * which returns a pointer to an inclusion anchor structure.
  1905.  * Call testing procedure with a pointer to this structure and test point
  1906.  * _point_, returns 1 if inside, 0 if outside.
  1907.  * Call cleanup with pointer to inclusion anchor structure to free space.
  1908.  *
  1909.  * CONVEX must be defined for this test; it is not usable for general polygons.
  1910.  */
  1911.  
  1912. #ifdef    CONVEX
  1913. /* make inclusion wedge set */
  1914. pInclusionAnchor InclusionSetup( pgon, numverts )
  1915. double    pgon[][2] ;
  1916. int    numverts ;
  1917. {
  1918. int    pc, p1, p2, flip_edge ;
  1919. double    ax,ay, qx,qy, wx,wy, len ;
  1920. pInclusionAnchor pia ;
  1921. pInclusionSet    pis ;
  1922.  
  1923.     /* double the first edge to avoid needing modulo during test search */
  1924.     pia = (pInclusionAnchor)malloc( sizeof( InclusionAnchor )) ;
  1925.     MALLOC_CHECK( pia ) ;
  1926.     pis = pia->pis =
  1927.         (pInclusionSet)malloc( (numverts+1) * sizeof( InclusionSet )) ;
  1928.     MALLOC_CHECK( pis ) ;
  1929.  
  1930.     pia->hi_start = numverts - 1 ;
  1931.  
  1932.     /* get average point to make wedges from */
  1933.     qx = qy = 0.0 ;
  1934.     for ( p2 = 0 ; p2 < numverts ; p2++ ) {
  1935.     qx += pgon[p2][X] ;
  1936.     qy += pgon[p2][Y] ;
  1937.     }
  1938.     pia->qx = qx /= (double)numverts ;
  1939.     pia->qy = qy /= (double)numverts ;
  1940.  
  1941.     /* take cross product of vertex to find handedness */
  1942.     pia->flip_edge = flip_edge =
  1943.         (pgon[0][X] - pgon[1][X]) * (pgon[1][Y] - pgon[2][Y] ) >
  1944.         (pgon[0][Y] - pgon[1][Y]) * (pgon[1][X] - pgon[2][X] ) ;
  1945.  
  1946.  
  1947.     ax = pgon[0][X] - qx ;
  1948.     ay = pgon[0][Y] - qy ;
  1949.     len = sqrt( ax * ax + ay * ay ) ;
  1950.     if ( len == 0.0 ) {
  1951.     fprintf( stderr, "sorry, polygon for inclusion test is defective\n" ) ;
  1952.     exit(1) ;
  1953.     }
  1954.     pia->ax = ax /= len ;
  1955.     pia->ay = ay /= len ;
  1956.  
  1957.     /* loop through edges, and double last edge */
  1958.     for ( pc = p1 = 0, p2 = 1
  1959.     ; pc <= numverts
  1960.     ; pc++, p1 = p2, p2 = (++p2)%numverts, pis++ ) {
  1961.  
  1962.     /* wedge border */
  1963.     wx = pgon[p1][X] - qx ;
  1964.     wy = pgon[p1][Y] - qy ;
  1965.     len = sqrt( wx * wx + wy * wy ) ;
  1966.     wx /= len ;
  1967.     wy /= len ;
  1968.  
  1969.     /* cosine of angle from anchor border to wedge border */
  1970.     pis->dot = ax * wx + ay * wy ;
  1971.     /* sign from cross product */
  1972.     if ( ( ax * wy > ay * wx ) == flip_edge ) {
  1973.         pis->dot = -2.0 - pis->dot ;
  1974.     }
  1975.  
  1976.     /* edge */
  1977.     pis->ex = pgon[p1][Y] - pgon[p2][Y] ;
  1978.     pis->ey = pgon[p2][X] - pgon[p1][X] ;
  1979.     pis->ec = pis->ex * pgon[p1][X] + pis->ey * pgon[p1][Y] ;
  1980.  
  1981.     /* check sense and reverse plane eqns if need be */
  1982.     if ( flip_edge ) {
  1983.         pis->ex = -pis->ex ;
  1984.         pis->ey = -pis->ey ;
  1985.         pis->ec = -pis->ec ;
  1986.     }
  1987.     }
  1988.     /* set first angle a little > 1.0 and last < -3.0 just to be safe. */
  1989.     pia->pis[0].dot = -3.001 ;
  1990.     pia->pis[numverts].dot = 1.001 ;
  1991.  
  1992.     return( pia ) ;
  1993. }
  1994.  
  1995. /* Find wedge point is in by binary search, then test wedge */
  1996. int InclusionTest( pia, point )
  1997. pInclusionAnchor    pia ;
  1998. double    point[2] ;
  1999. {
  2000. register double tx, ty, len, dot ;
  2001. int    inside_flag, lo, hi, ind ;
  2002. pInclusionSet    pis ;
  2003.  
  2004.     tx = point[X] - pia->qx ;
  2005.     ty = point[Y] - pia->qy ;
  2006.     len = sqrt( tx * tx + ty * ty ) ;
  2007.     /* check if point is exactly at anchor point (which is inside polygon) */
  2008.     if ( len == 0.0 ) return( 1 ) ;
  2009.     tx /= len ;
  2010.     ty /= len ;
  2011.  
  2012.     /* get dot product for searching */
  2013.     dot = pia->ax * tx + pia->ay * ty ;
  2014.     if ( ( pia->ax * ty > pia->ay * tx ) == pia->flip_edge ) {
  2015.     dot = -2.0 - dot ;
  2016.     }
  2017.  
  2018.     /* binary search through angle list and find matching angle pair */
  2019.     lo = 0 ;
  2020.     hi = pia->hi_start ;
  2021.     while ( lo <= hi ) {
  2022.     ind = (lo+hi)/ 2 ;
  2023.     if ( dot < pia->pis[ind].dot ) {
  2024.         hi = ind - 1 ;
  2025.     } else if ( dot > pia->pis[ind+1].dot ) {
  2026.         lo = ind + 1 ;
  2027.     } else {
  2028.         goto Foundit ;
  2029.     }
  2030.     }
  2031.     /* should never reach here, but just in case... */
  2032.     fprintf( stderr,
  2033.         "Hmmm, something weird happened - bad dot product %lg\n", dot);
  2034.  
  2035.     Foundit:
  2036.  
  2037.     /* test if the point is outside the wedge's exterior edge */
  2038.     pis = &pia->pis[ind] ;
  2039.     inside_flag = ( pis->ex * point[X] + pis->ey * point[Y] <= pis->ec ) ;
  2040.  
  2041.     return( inside_flag ) ;
  2042. }
  2043.  
  2044. void InclusionCleanup( p_inc_anchor )
  2045. pInclusionAnchor p_inc_anchor ;
  2046. {
  2047.     free( p_inc_anchor->pis ) ;
  2048.     free( p_inc_anchor ) ;
  2049. }
  2050. #endif
  2051.  
  2052.  
  2053. /* ======= Crossings Multiply algorithm =================================== */
  2054.  
  2055. /*
  2056.  * This version is usually somewhat faster than the original published in
  2057.  * Graphics Gems IV; by turning the division for testing the X axis crossing
  2058.  * into a tricky multiplication test this part of the test became faster,
  2059.  * which had the additional effect of making the test for "both to left or
  2060.  * both to right" a bit slower for triangles than simply computing the
  2061.  * intersection each time.  The main increase is in triangle testing speed,
  2062.  * which was about 15% faster; all other polygon complexities were pretty much
  2063.  * the same as before.  On machines where division is very expensive (not the
  2064.  * case on the HP 9000 series on which I tested) this test should be much
  2065.  * faster overall than the old code.  Your mileage may (in fact, will) vary,
  2066.  * depending on the machine and the test data, but in general I believe this
  2067.  * code is both shorter and faster.  This test was inspired by unpublished
  2068.  * Graphics Gems submitted by Joseph Samosky and Mark Haigh-Hutchinson.
  2069.  * Related work by Samosky is in:
  2070.  *
  2071.  * Samosky, Joseph, "SectionView: A system for interactively specifying and
  2072.  * visualizing sections through three-dimensional medical image data",
  2073.  * M.S. Thesis, Department of Electrical Engineering and Computer Science,
  2074.  * Massachusetts Institute of Technology, 1993.
  2075.  *
  2076.  */
  2077.  
  2078. /* Shoot a test ray along +X axis.  The strategy is to compare vertex Y values
  2079.  * to the testing point's Y and quickly discard edges which are entirely to one
  2080.  * side of the test ray.  Note that CONVEX and WINDING code can be added as
  2081.  * for the CrossingsTest() code; it is left out here for clarity.
  2082.  *
  2083.  * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
  2084.  * _point_, returns 1 if inside, 0 if outside.
  2085.  */
  2086. int CrossingsMultiplyTest( pgon, numverts, point )
  2087. double    pgon[][2] ;
  2088. int    numverts ;
  2089. double    point[2] ;
  2090. {
  2091. register int    j, yflag0, yflag1, inside_flag ;
  2092. register double    ty, tx, *vtx0, *vtx1 ;
  2093.  
  2094.     tx = point[X] ;
  2095.     ty = point[Y] ;
  2096.  
  2097.     vtx0 = pgon[numverts-1] ;
  2098.     /* get test bit for above/below X axis */
  2099.     yflag0 = ( vtx0[Y] >= ty ) ;
  2100.     vtx1 = pgon[0] ;
  2101.  
  2102.     inside_flag = 0 ;
  2103.     for ( j = numverts+1 ; --j ; ) {
  2104.  
  2105.     yflag1 = ( vtx1[Y] >= ty ) ;
  2106.     /* Check if endpoints straddle (are on opposite sides) of X axis
  2107.      * (i.e. the Y's differ); if so, +X ray could intersect this edge.
  2108.      * The old test also checked whether the endpoints are both to the
  2109.      * right or to the left of the test point.  However, given the faster
  2110.      * intersection point computation used below, this test was found to
  2111.      * be a break-even proposition for most polygons and a loser for
  2112.      * triangles (where 50% or more of the edges which survive this test
  2113.      * will cross quadrants and so have to have the X intersection computed
  2114.      * anyway).  I credit Joseph Samosky with inspiring me to try dropping
  2115.      * the "both left or both right" part of my code.
  2116.      */
  2117.     if ( yflag0 != yflag1 ) {
  2118.         /* Check intersection of pgon segment with +X ray.
  2119.          * Note if >= point's X; if so, the ray hits it.
  2120.          * The division operation is avoided for the ">=" test by checking
  2121.          * the sign of the first vertex wrto the test point; idea inspired
  2122.          * by Joseph Samosky's and Mark Haigh-Hutchinson's different
  2123.          * polygon inclusion tests.
  2124.          */
  2125.         if ( ((vtx1[Y]-ty) * (vtx0[X]-vtx1[X]) >=
  2126.             (vtx1[X]-tx) * (vtx0[Y]-vtx1[Y])) == yflag1 ) {
  2127.         inside_flag = !inside_flag ;
  2128.         }
  2129.     }
  2130.  
  2131.     /* Move to the next pair of vertices, retaining info as possible. */
  2132.     yflag0 = yflag1 ;
  2133.     vtx0 = vtx1 ;
  2134.     vtx1 += 2 ;
  2135.     }
  2136.  
  2137.     return( inside_flag ) ;
  2138. }
  2139.