Rcpp point in polygon

gistfile1.cpp

#include <Rcpp.h>

using namespace Rcpp;
// ref: http://geomalgorithms.com/a03-_inclusion.html
// Copyright 2000 softSurfer, 2012 Dan Sunday
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.
 

// a Point is defined by its coordinates {int x, y;}
//===================================================================
 

// isLeft(): tests if a point is Left|On|Right of an infinite line.
//    Input:  three points P0, P1, and P2
//    Return: >0 for P2 left of the line through P0 and P1
//            =0 for P2  on the line
//            <0 for P2  right of the line
//    See: Algorithm 1 "Area of Triangles and Polygons"
  int isLeft( NumericVector P0, NumericVector P1, NumericVector P2 )
	{
	    return ( (P1[0] - P0[0]) * (P2[1] - P0[1])
	            - (P2[0] -  P0[0]) * (P1[1] - P0[1]) );
	}
//===================================================================


// cn_PnPoly(): crossing number test for a point in a polygon
//      Input:   P = a point,
//               V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
//      Return:  0 = outside, 1 = inside
// This code is patterned after [Franklin, 2000]

int	cn_PnPoly( NumericVector P, NumericMatrix V) // P is vector of 2, (x, y). V is nx2 matrix. 
	{
		int n = V.nrow() - 1; 
	    int    cn = 0;    // the  crossing number counter

	    // loop through all edges of the polygon
	    for (int i=0; i<n; i++) {    // edge from V(i]  to V(i+1]
	       if (((V(i, 1) <= P[1]) && (V(i+1, 1) > P[1]))     // an upward crossing
	        || ((V(i, 1) > P[1]) && (V(i+1, 1) <=  P[1]))) { // a downward crossing
	            // compute  the actual edge-ray intersect x-coordinate
	            float vt = (float)(P[1]  - V(i,1)) / (V(i+1, 1) - V(i, 1));
	            if (P[0] <  V(i, 0) + vt * (V(i+1, 0) - V(i, 0))) // P[0] < intersect
	                 ++cn;   // a valid crossing of y=P[1] right of P[0]
	        }
	    }
	    return (cn % 1);    // 0 if even (out), and 1 if  odd (in)
	}
//===================================================================


// wn_PnPoly(): winding number test for a point in a polygon
//      Input:   P = a point,
//               V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
//      Return:  wn = the winding number (=0 only when P is outside)
// [[Rcpp::export]]
int wn_PnPoly( NumericVector P, NumericMatrix V)
	{
		int n = V.nrow()-1, wn=0;  // the  winding number counter

	    // loop through all edges of the polygon
	    for (int i=0; i<n; i++) {   // edge from V[i] to  V[i+1]
	        if (V(i, 1) <= P[1]) {          // start y <= P.y
	            if (V(i+1, 1)  > P[1])      // an upward crossing
	                 if (isLeft( V(i, _) , V(i+1, _), P) > 0)  // P left of  edge
	                     ++wn;            // have  a valid up intersect
	        }
	        else {                        // start y > P.y (no test needed)
	            if (V(i+1, 1)  <= P[1])     // a downward crossing
	                 if (isLeft( V(i, _), V(i+1, _), P) < 0)  // P right of  edge
	                     --wn;            // have  a valid down intersect
	        }
	    }
	    return wn;
	}
//===================================================================

//' @Reference \url{http://geomalgorithms.com/a03-_inclusion.html}
//' @return The point is outside/on board only when this "winding number" wn = 0; otherwise, the point is inside.
// [[Rcpp::export]]
IntegerVector PnPoly( NumericMatrix P, NumericMatrix V)
{
  int n = P.nrow();
  IntegerVector pos(n);
  for(int i =0; i < n; i++)
  {
    pos[n] = wn_PnPoly(P(i,_), V);
  }
  return pos; 
}
asf 
//' @rdname GIS
//' @reference \url{http://geomalgorithms.com/a13-_intersect-4.html#intersect2D_SegPoly()}