Point in polygon

Central Almaty has 7909 crossings, and the city has eight districts. Each crossing needs its district: who repairs the traffic light, who clears the snow. The program has only the coordinates of a crossing and the district boundaries, polylines of hundreds of points. Below all the crossings are coloured by district, and every colour was set by the same test: a ray from the point to the east and a count of how many times it crossed the boundary. Tap the map anywhere.

The test assigned 7445 crossings to districts, and not one landed in two districts at once. Another 459 stand exactly on a boundary, to the metre, and five landed in no district. Neither is a fault of the test. District boundaries in OpenStreetMap are in many places drawn right through crossings, and on wide avenues neighbouring districts split the street not along one line, leaving a no-man’s strip between their boundaries. All this is below; let us start with why the ray works at all.

The idea

Walk from the point in a straight line in one direction until you are out of town. Each crossing of the district boundary switches “inside” to “outside” and back, and at the end of the walk you are certainly outside. So if there was an odd number of crossings, you started inside.

Draw a ray from the point and count how many times it crosses the boundary of the polygon. An odd number means the point is inside, an even one outside. Side AB crosses the eastward ray from point P if exactly one of its ends lies above P and the crossing point is east of P; the latter is checked by the sign of a cross product, without division.

This is the even-odd rule, or ray test. The direction of the ray is arbitrary; east is most convenient: the crossing with a side is sought on the line y = P_y. A side from A to B crosses it if one end is above and the other is not. And whether the crossing is east of the point is told by the sign of the cross product cross(A, B, P): if the side goes up, the crossing is to the east when P is to its left; if down, when to its right. No fractions, no division: for integer coordinates the answer is exact.

Why it is correct

Behind the test is the Jordan curve theorem (1887): a simple closed curve divides the plane into two parts, the inside and the outside. In general its proof is hard, but for a polygon the ray itself provides it. Move the point without touching the boundary. The number of crossings of its ray changes only when the ray passes through a vertex, and then it changes by 0 or by 2: the two sides at the vertex either replace one another or appear on the ray and leave it together. So the parity is constant on any piece of the plane that can be traversed without crossing the boundary. When the point crosses a side, its ray loses or gains exactly one crossing, and the parity changes. So the points fall into two classes, and far from the polygon there are zero crossings: that is “outside”. Oswald Veblen in 1905 reproached Jordan precisely for this: he had accepted the theorem for polygons without proof.

A ray through a vertex

The first trap is a ray that passes exactly through a vertex. If a side is counted when the point’s y lies in the segment [min, max] of its ends inclusive, the ray touches both sides at the vertex. When the boundary passes through the ray right across, that is one crossing, but two are counted, and the answer flips. On a map of city districts this is not a rarity but the rule: the coordinates are integers, the boundaries are drawn through crossings, and rays from neighbouring crossings keep running into vertices. Switch the scene to “Naive ray”: such counting errs on 355 crossings of 7909, every twenty-second one.

The cure is one condition: a side is counted if exactly one of its ends is strictly above the ray. A vertex on the ray then counts as if it were slightly below it. If the boundary passes through the ray right across, one side at the vertex goes up and is counted, the other goes down and is not: one crossing. If the boundary only touches the ray, both sides are on one side of it: both are counted or neither, and the parity does not change. A horizontal side on the ray is never counted: its neighbours take care of themselves. Try it on the crown below: point P stands at the height of the dips between the teeth, and the right flank of the crown bends exactly at this height.

A point on the boundary

The second trap is a point on a side itself. The ray from it starts on the boundary, and the parity means nothing: the answer depends on how the rule rounds the crossing at the point itself. So “on the boundary” is checked separately and first of all: P lies on side AB if cross(A, B, P) = 0 and P is inside the rectangle spanned by A and B. Almaty has 459 such crossings: 222 lie on the boundary of one district, 235 on the shared boundary of two, and two in a corner where three districts meet. All 459 are vertices of boundaries: in OpenStreetMap the district boundary was in many places drawn right through road nodes.

If neighbouring polygons share sides exactly, points on the boundary can be got rid of by convention. Graphics cards do this: a pixel whose centre lies on the shared side of two triangles is drawn by exactly one of them; in Direct3D this is the top-left rule. With Almaty’s districts that will not work: neighbouring districts split wide avenues not along one line. At the intersection of Abay and Abylai Khan avenues the boundary of Almaly district passes through the crossing with the northern carriageway of Abay, and the boundary of Bostandyk district through the crossing with the southern one, fifteen metres away. The median strip between them belongs to no district. The five crossings the test gave to nobody stand on such strips; choose the example “A no-man’s strip on Abay Avenue” in the scene. The test is right here; it is the data that is imprecise.

The winding number

The ray has a second way of counting. Let us distinguish the direction in which the boundary crosses the ray: from below upwards, +1; from above downwards, −1. The sum is the winding number of the boundary around the point: how many times it goes around the point anticlockwise minus the times clockwise. If you stand at the point and keep looking at someone running along the boundary, the direction of your gaze passes through east exactly when the runner crosses the ray, and turns anticlockwise if they are running upwards. Every full turn passes through east once.

For a simple polygon the winding number inside is +1 or −1, depending on the direction of traversal, and outside 0, and the answer is the same as with parity: at all of Almaty’s crossings, except those on a boundary, both tests agree. The difference shows when the boundary crosses itself. Choose the pentagram in the scene: from its centre the ray crosses two sides, an even number, “outside”, but crosses both from below upwards, winding number 2, “inside”. Hence the two fill rules found in SVG and on the browser canvas: evenodd fills a point with an odd number of crossings, nonzero one with a non-zero winding number. The default in both is nonzero. The loop in the same scene goes around its middle twice, and around the right ear clockwise: drag the point there, and the winding number becomes −1.

The winding number used to be computed as the sum of the angles at which the sides are seen from the point, through arctangents. Counting crossings with signs gives the same integer without trigonometry. This is the variant analysed by Dan Sunday (2001) and by Kai Hormann and Alexander Agathos in a paper of the same year on the point in a polygon with self-intersections.

Code

Three tests on small shapes: the ray with the half-open rule and a boundary check, the naive count, and the winding number. The crown is the same as in the scene: the naive count errs for three points at height 3, and for a point in a dip answers “inside”, although it is on the boundary. The pentagram is for comparing parity and winding.

And this is the first scene in numbers: all crossings against all districts, counting how many sides had to be checked.

Many points: strips

The ray test looks at all n sides, although only those spanning P_y in height can cross the ray. If there are many points, the sides are worth laying out in advance. Cut the district into horizontal strips of equal height, as many strips as sides, and record each side in all the strips it touches. The ray from a point lies in one strip, and only its sides need looking at: any side that can cross the ray is there. In Almaty’s districts a side is tens of metres of boundary, and a strip usually holds a handful. The listing above counted 16,429,057 side checks without strips and 61,049 with them, 269 times fewer. Preparation is one pass over the sides.

Trying all sides grows linearly with the number of sides, strips far more slowly: at ten thousand sides they are hundreds of times faster in Python and dozens in Go. The weak spot is the shape. If many sides span the whole polygon in height, like the teeth of a comb, they fall into every strip, and the gain disappears. A guarantee of O(log n) per point for any shape is given by structures for planar subdivisions: the slabs of Dobkin and Lipton (1976), which cut the plane through every vertex and search by binary search twice, and their memory-saving descendants. And if the points come on their own, without a polygon at hand, they are first gathered into a k-d tree or an R-tree and asked which of them fell into a district’s bounding box, and the exact test is done only for those.

What it costs

One point is O(n) time and O(1) memory: one check per side, with two differences and two multiplications in it. For a simple polygon it cannot be done faster without preparation: without looking at some side, you cannot know whether it passes between the point and the rest of the boundary. With strips it is O(n + S) to prepare, where S is how many times sides went into strips, and O(1 + k) per point, where k is the sides in its strip.

Precision is as for all geometry on cross products. In integers the answer is exact: with coordinates up to 10⁹ the differences are up to 2·10⁹ and the products up to 4·10¹⁸, and everything fits in int64. That is why the crossings and boundaries in this article are converted to integer metres. In floating point, points right at the boundary may end up on either side of it: the same trap as with turning left or right. If the answer “on the boundary” matters, coordinates are converted to integers or computed with a tolerance, honestly answering “too close to the boundary”.

Where it works

• Maps and delivery. Which district a house is in, whether an address falls into a delivery zone, which cadastral plot is under a point: the ray test after a rough filter by bounding box, as in the listing.
• Vector graphics. The fill-rule fill rules in SVG and the canvas fill() are the parity of crossings and the winding number. The letter “O” is two contours, and its hole stays a hole by either rule if the inner contour is traversed in the opposite direction: the winding number there is 1 − 1 = 0, and there are two crossings.
• Interfaces and games. Whether a click hit a country on a map or a shape on the screen. In this article a district is selected this way in both scenes.
• Rasterisation. A graphics card decides “is the pixel centre inside” for each triangle, and the top-left rule what to do with pixels on a shared side.

History

In August 1962 the journal Communications of the ACM printed Moshe Shimrat’s procedure “Algorithm 112: Position of point relative to polygon”, the earliest publication of the ray test usually cited. In December the journal published a certification of it by another author: at the time procedures from the algorithms section were customarily rechecked. Since then the test has been rediscovered many times, and people stumbled on the same vertices. In 1990 the Ray Tracing News mailing list described a solution that, when the ray hit a vertex, proposed simply firing another ray. Mark VandeWettering replied that it is simpler to agree to treat a vertex on the ray as lying on one side of it, and the answer comes out right. Asked whose idea it was, he wrote that Eric Haines had shown it to him, and Haines did not claim authorship either. The winding number without trigonometry was described by Dan Sunday in 2001; the same year Kai Hormann and Alexander Agathos analysed the point in a polygon of arbitrary shape, with self-intersections, in the journal Computational Geometry.

Problems

In the first problem crossings have to be sorted into “inside”, “outside” and “on the boundary”, with vertices on the ray, horizontal sides and coordinates up to a billion; two tests are real districts of Almaty. In the second, count how many times a runner’s track went around each point: parity is not enough there, the winding number is needed.