4
$\begingroup$

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?

$\endgroup$
4
$\begingroup$

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).

$\endgroup$
  • $\begingroup$ I was, indeed, referring to the cases where the ray crosses vertices or lie collinear to an edge. I agree, that the above algorithm could have "gotchas" were the point near an edge or vertex. In my implementation, all of which I do not show above, the function includes a distance_to_edge tolerance $s_toler$. Should a point find itself with distance from the edge within this value the point is instead declared to lie on the boundary. As long as this tolerance $s_toler > 0$ near-edge artifacts are largely ignored. $\endgroup$ – Mountainview Mar 24 '16 at 19:10
  • $\begingroup$ There's a few things I omitted in my answer. 1. If you have doubts about the stability of the algorithm, you can cast multiple rays to increase confidence in the result. 2. I mentioned using quadtrees, but using that makes the theoretical bound jump to O(n log(n)). 3. More on the software side, the line integral uses atan2. General intuition is that antan2 would be slower. Maybe consider using complex numbers multiplication to sum angles and count crossings? 4. Also software, profiling helps here. Are there convex shapes? Self intersecting? How many lines does a ray cross in average? $\endgroup$ – ZeroUltimax Mar 24 '16 at 19:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.