physics: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.

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*, *P _{y}* is the y-coordinate and

$$ %\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} $$

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 *P _{y}* 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 *P _{x}* 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)