Segment intersection

Two streets on a map cross. If it is a junction, in OpenStreetMap they have a shared node: a car can turn from one onto the other. If there is no shared node, one street passes over the other, on a bridge or in a tunnel, and such a road must carry a layer tag. Below are all the streets of central Almaty: 6275 lines, 37,316 segments between neighbouring points. The program found all the places where two segments meet without a shared node. Tap a circle.

There are almost 700 million pairs of segments here, and the program looked at only 146,150 of them, the pairs that landed in the same grid cell, and checked exactly 52,018 of those, whose bounding boxes overlap. There were 166 crossings without a shared node, and every single one is a bridge or a tunnel: the crossing streets are on different layers. There is not a single place in the city centre where streets cross on one layer and forgot a shared node. But two oddities of another kind turned up: a piece of road drawn twice, two OpenStreetMap ways on the same two nodes, and two different nodes of Ryskulbekov Street a centimetre and a half apart. To find all this it is not enough to know whether two segments cross: you need to understand exactly how.

The idea

Whether segments AB and CD cross is decided by four cross products: the ends of each segment must lie on different sides of the other. This is covered in the article on vectors and orientation. Here is the next step: where exactly they meet, and what to do when there are thousands of segments.

Write the segments as A + t·r and C + u·s, where r = B − A, s = D − C, and t and u run from 0 to 1. Taking the cross product of the equation A + t·r = C + u·s with s gives t = (C − A) × s / (r × s), and with r, u = (C − A) × r / (r × s). The segments meet if both numbers are between 0 and 1. The denominator r × s is zero when the segments are parallel, which is a separate case.

For integer coordinates t and u are fractions with a common integer denominator, and they can be compared with 0 and 1 without dividing: 0 ≤ t_num ≤ denominator. The crossing point A + t·r is a fraction too. Where t and u lie also shows how the segments meet: both strictly inside, a cross; one at an end, a letter T; both at ends, a touch end to end.

All the cases

Drag the ends and check each case. The crossing point under the picture is written exactly, in fractions, with a decimal approximation beside it.

When r × s = 0, the segments are parallel. If point C is not on line AB, they do not meet. If it is, all four points are on one line, and the question reduces to segments on a number line: order the ends along the line by x (or by y for a vertical one) and compare. A common piece of positive length is an overlap, one common point is a touch, otherwise a gap. An overlap is not a point but a segment, and a program that always returns “the crossing point” gets it wrong. In the Almaty data there is one such overlap: that very piece of road drawn twice.

The point is a fraction, and a big one

The check “do they meet” is exact in integers, but the crossing point is a fraction, and its numerator grows fast. With coordinates up to 10⁹ the denominator r × s reaches 8·10¹⁸, on the edge of int64, and the numerator of a point’s coordinate up to 10²⁸. In the listing below, the numerator of y for segments with coordinates around a billion is a 27-digit number. Python handles this by itself; Go needs math/big.

In practice the approach is this: decide “do they meet and how” exactly, in integers, and compute the point itself, if it is needed for drawing, in double. The other way round does not work: if you find the point in floating point and then ask whether it lies on both segments, at the edge the answer is decided by rounding error, the same trap as with turning left or right. Maps run into this all the time: among Almaty’s streets converted to decimetres, 47,508 pairs of segments touch at shared nodes, and every such touch is exactly at an end, at t = 0 or t = 1.

Code

The crossing point in fractions and all the cases from the scene, and at the end segments with coordinates around a billion.

And this is the first scene in numbers: all the streets of central Almaty, a grid of 50 × 50 m cells and an analysis of each meeting.

Many segments: a grid

Checking all pairs of n segments is n(n − 1)/2 checks: 696,223,270 for 37,316 segments. But only segments that lie close together can meet. Cut the map into square cells and record each segment in all the cells its bounding box touches. Two segments that have a common point necessarily land in at least one common cell, the one where that point lies. So only pairs within cells need checking.

One subtlety: a long segment lies in several cells, and a pair may meet in more than one of them. To avoid checking it twice, agree to check a pair only in the cell containing the lower-left corner of the intersection of their bounding boxes: there is only one such cell. Pairs with disjoint bounding boxes are discarded along the way. In the listing above 146,150 pairs ended up in shared cells, almost five thousand times fewer than all pairs, and 52,018 of them needed the exact check.

The cell size is a choice. Move the slider in the first scene. With small cells there are few pairs, but each segment is recorded in many cells, and the work goes into laying them out: with 10 m cells there are more than ten times as many cell entries as segments. With large cells there are few entries, but hundreds of segments in each: with a kilometre cell there are more than six million pairs in cells, forty times as many as with 50 m. The number of exact checks is always the same: the cell size decides only how much excess must be discarded. A good cell is of the order of a typical segment’s length. A grid is a hash table keyed by cell number, and it has the same weakness: if segments crowd into one cell or some are very long, there are no guarantees, and in the worst case it is O(n²) again.

The sweep line gives a guarantee. Michael Shamos and Dan Hoey showed in 1976 how to find out in O(n log n) whether there is at least one crossing pair among n segments, and Jon Bentley and Thomas Ottmann in 1979 how to find all k crossings in O((n + k) log n). The line moves from left to right, keeps the segments it currently crosses in top-to-bottom order, and checks only neighbours in this order. Only segments that are neighbours somewhere can cross. A separate article is devoted to this.

What it costs

One pair is O(1): a few multiplications and comparisons. n segments by brute force are O(n²) pairs. A grid is O(n + P), where P is the pairs in shared cells; on streets that is hundredths of a per cent of all pairs, but in the worst case P grows as n². The sweep line is O((n + k) log n); the factor k cannot be removed at all, since all crossings must at least be listed, and there can be up to n²/2 of them. The optimal time O(n log n + k) is achieved by the algorithm of Bernard Chazelle and Herbert Edelsbrunner (1992), but it is complex, and in practice a grid or a sweep line is enough.

Where it works

• Map validation. The JOSM editor’s validator looks for roads, rivers and railways that cross on the same layer without a shared node, and suggests setting a layer, and perhaps a bridge or tunnel. That is exactly what the first scene does; no such errors were found in central Almaty.
• Polygons. Before computing an area or checking a point inside, one checks that the sides of the contour do not cross, otherwise both formulas answer the wrong question.
• Games and robots. Did a bullet hit a wall, does a guard see the player, will a robot hit an obstacle on its way: the crossing of a path segment with wall segments, usually with the same grid to select candidates.
• Drawings and boards. Traces on a printed circuit board must not cross, and a part’s contours must not self-intersect: programs check this before manufacturing.

Problems

In the first problem you must find where two segments meet, exactly, in fractions, with overlaps and coordinates up to a billion, when the numerator no longer fits in 64 bits. In the second, count the bridges: how many pairs of segments among tens of thousands cross without a shared node; trying all pairs will not fit in the time.