Point In Polygon

/*
 * Algorithm : Point in Polygon By ray shooting
 *             Adapted from "Computational Geometry in C" (Second Edition)
 * Complexity : O( N )
 * Note : For each query point Pt, InPoly returns one of four char's:
 *  i : Pt is strictly interior to Pl
 *  o : Pt is strictly exterior to Pl
 *  v : Pt is a vertex of Pl
 *  e : Pt lies on the relative interior of an edge of Pl
 *
*/
#include<stdio.h>
#include<math.h>
#define X 0
#define Y 1
#define DIM 2
typedef long POINT[DIM];   /* type integer point */
char PointInPoly( POINT Pt,POINT *Pl, long N ){
    long i,j,d;
    long Rcross = 0; /* number of right edge/ray crossings */
    long Lcross = 0; /* number of left edge/ray crossings */

    /* Shift so that Pt is the origin. Note this destroys the polygon.*/
    for( i=0;i<N;i++ ){
        for( d=0;d<DIM;d++ ) Pl[i][d] = Pl[i][d] - Pt[d];
    }
    /* For each edge e=(i-1,i), see if crosses ray. */
    for( i=0;i<N;i++ ){
        /* First see if Pt=(0,0) is a vertex. */
        if( Pl[i][X]==0 and Pl[i][Y]==0 ) return 'v';
        j = (i+N-1)%N;

        bool Rstrad = (Pl[i][Y]>0) != (Pl[j][Y]>0);
        bool Lstrad = (Pl[i][Y]<0) != (Pl[j][Y]<0);

        /* if e "straddles" the x-axis... */
        if( Rstrad or Lstrad ){
            /* e straddles ray, so compute intersection with ray. */
            double x = ( Pl[i][X]*(double)Pl[j][Y] - Pl[j][X]*(double)Pl[i][Y]) / (double)(Pl[j][Y]-Pl[i][Y]);
            /* crosses ray if strictly positive intersection. */
            if( Rstrad and x>0 ) Rcross++;
            if( Lstrad and x<0 ) Lcross++;
        }
    }
    /* Pt on the edge if left and right cross are not the same parity. */
    if( (Rcross%2) != (Lcross%2) ) return 'e';
    /* Pt inside iff an odd number of crossings. */
    if( (Rcross%2) == 1 ) return 'i';
    else return 'o';
}