User Tools

Site Tools


Point-in-Polygon Problem with Simulation of Simplicity

This article discusses the Simulation of Simplicity (SoS) technique to find a simple, fast and robust algorithm to solve the Point-in-Polygon problem.

One of the most effective ways to determine, if a point lies inside or outside a given polygon is to use the Jordan curve theorem. In this case, a ray is drawn from the inspected point P to infinity. If this ray intersects the polygon odd number times, then the point lies inside it. If the number of intersections is even, the point lies outside the polygon. Without affecting generality, one can limit themselves only to horizontal rays pointing to positive infinity.

In this case, the Jordan curve is a polygon. If the testing ray intersects this polygon in a vertex or multiple vertices, the result may be wrong, when a computer is used. It is caused by the fact that there is only a finite set of numbers in computers. The probability of hitting a vertex is higher than zero and depends on the chosen accuracy. When a vertex or vertices are hit, it is called a degenerate case. Cases marked B, D, E, and F in the following figure are degenerate.

Point-in-polygon problem - degenerate cases
Fig 1: Point-in-polygon problem. Cases A and C are non-degenerate, B, D, E, and F are degenerate cases.

Since the degenerations may lead to incorrect results, one approach is to check for these degenerations. If the ray hits a vertex, draw another one. This may result in bugs, if one or more degenerate cases are not covered, and the checks cost some machine time. If the point-in-polygon detection is run many times, the slow-down due to the checks may be significant.

Another approach is to use the Simulation of Simplicity. It is a technique that introduces an infinitesimal virtual perturbation, which does not affect the result, but eliminates the degenerations.

The test, whether the ray intersects the edges of a polygon can be divided into a set of tests, whether this ray intersects an individual edge. This is in turn repeated for all edges. The standard way to determine this is to parametrize both - the ray and the edge. For the horizontal ray, the following parametrization can be written:

$$ \eqalign{ x(t) &= P_x + t,\cr y(t) &= P_y\tag{1}, } $$

where $P_x$ is the x-coordinate of the inspected point P, Py is the y-coordinate and t is the parameter. The edge can be parametrized analogously to (1):

$$ %\displaystyle \eqalign{ x(u) &= (V_{2x} - V_{1x})u + V_{1x},\cr y(u) &= (V_{2y} - V_{1y})u + V_{1y},\tag{2} } $$

where $V_{1x}$ … $V_{2y}$ are the coordinates of the corresponding vertices. The ray intersects with the edge, if all the following conditions are met:

$$ \displaystyle \eqalign{ x(t) &= x(u), \cr y(t) &= y(u), \cr t &> 0, \cr 0 \leq &u \leq 1. }\tag{3} $$

Intersection of the testing ray and a polygon edge
Fig 2: Intersection of the testing ray and a polygon edge.

Fig 2 shows that if the ray and the edge intersect, the inspected point P must lie within the filled area (band). The condition for Py thus is:

$$ \eqalign{ V_{2y} \geq &P_y \geq V_{1y}, {\quad\rm or}\cr V_{1y} \geq &P_y \geq V_{2y}. }\tag{4} $$

However, if $P_y = V_{1y}$ or $P_y = V_{2y}$, the result may be wrong, since this is a degeneration. To avoid this, one can shift the edge (both vertices) vertically by an infinitesimal distance ε so that:

$$ \eqalign{ \varepsilon &> 0,\cr \varepsilon &\to 0. }\tag{5} $$

The condition (4) now is:

$$ \displaystyle\eqalign{ V_{2y}+\varepsilon \geq &P_y \geq V_{1y}+\varepsilon, {\quad\rm or}\cr V_{1y}+\varepsilon \geq &P_y \geq V_{2y}+\varepsilon. }\tag{6} $$

Since ε is infinitesimally small, it can be vanished and the conditions are:

$$ \displaystyle\eqalign{ V_{2y} \geq &P_y > V_{1y}, {\quad\rm or}\cr V_{1y} \geq &P_y > V_{2y}. }\tag{7} $$

There are no more degenerate cases, no checks are needed.

From the first equation in (3), and the set of equations (1) and (2), the condition for Px may be calculated:

$$ \displaystyle P_x < \frac{(V_{2x}-V_{1x})(P_y-V_{1y})}{V_{2y}-V_{1y}}+V_{1x}.\tag{8} $$

There are many implementations of this test. The following one in C is used in ARTIMAGEN:

  int i, j, c = 0;
  for (i = 0, j = number_of_vertices-1; i < number_of_vertices; j = i++) {
    if ( ((vertices[i].y>p.y) != (vertices[j].y>p.y)) &&
         (p.x < (vertices[j].x-vertices[i].x) * (p.y-vertices[i].y) / (vertices[j].y-vertices[i].y) + vertices[i].x) )
       c = !c;
  return c;

p.x and p.y are coordinates of the inspected point. This code is pretty fast. It significantly helped to speed-up the ARTIMAGEN library.

For details, see the original article on arXiv.

physics/point-in-polygon_problem_with_simulation_of_simplicity.txt · Last modified: 2017/05/16 11:10 (external edit)