Some twenty years prior I was given the task to solve the Point in Polygon problem for a piece of commercial software. I solved it invoking the ray casting algorithm. After a variety of enhancements the function worked well enough, but I felt there was a better way. I found inspiration in Gauss's Law. Charges and the electric field aside, it boils down to, in 2D, a (contour) line integral. In a plane, the integral $\int_C d\theta = 2\pi$ if the point lies within the contour, and is zero if the point lies outside. With this in mind I reformulated the algorithm for arbitrary polygons:
double sum_theta = 0`
for( int i = 0 ; i < n_vtx ; i++ )
{
int j = (i + 1) % n_vtx;
XYPT p_i = ngon[i];
XYPT p_j = ngon[j];
XYPT vector_i = p_i - p_o; // vector from point o to vertex i
XYPT vector_j = p_j - p_o; // vector from point o to vertex j
double cp = vector_i.cross( vector_j ) // cross product
double dp = vector_i.dot ( vector_j ) // dot product
sum_theta+= Atan2( cp, dp );
}
if( sum > pi ) return INSIDE;
else return OUTSIDE;
I'm not asking if this works. It does, rather well. (And I see no reason to replace it; I like its simplicity and elegance.) My question to this forum is largely of curiosity: Indeed, how much slower is this method than doing it by ray-casting (which can present gotchas) and intersection count; or by counting winding number (WN)? The application this services never has self-intersecting polygons (so, it seems, WN has no particular advantage) and the vertex count is never excessive (n_max ~ 25).
Any thoughts?