Point in Polygon Problem: Has anybody invoked the line integral?

Question

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?

Answer 1

Summing the angles around a point is equivalent to calculating the winding number around that point. It will give you a correct answer in the general case.

Going more in depth, I'll try to show that your proposed algorithm has many drawbacks.

Algorithmic runtime

Algorithmically, your proposed algorithm must iterate over all $N$ points. Because of that, we can assume it runs in $\Theta(N)$ time.

In contrast, a ray casting algorithm checks against edges, of which there are $N$. The algorithm can be optimized using grids or quadtrees to reduce the amount of calculations needed. Since the ray might intersect all edges, that algorithm runs in $\mathcal{O}(N)$ time.

We can conclude that the ray casting algorithm will never perform worse than the angle sum.

Gotchas

You mentioned the ray casting algorithm as having gotchas. I'll assume you refer to the problems occurring when the ray crosses a vertex or an edge directly. First, your proposed algorithm has the same problem. What if the point you choose is on a vertex or an edge?

The solution to the problem for both algorithms is to define strict rules. That is, from the start, decide (or parametrize the algorithm) so that points on edges or vertices of the polygon are either all included or all excluded.

As an aside, the solution for ray casting's problem of intersecting vertices and edges can be solved using a similar rule. When a vertex is on the ray, it must be considered to belong to one side of the ray (more here).