Vectors and orientation

The navigator says: “In two hundred metres, turn right.” How does it know it is right? It has no eyes, only coordinates: a route is a polyline of points. Below is a real route through the streets of Almaty, from Republic Square to the Green Bazaar, the same as in the article on Dijkstra. The program got all the prompts on the right from the coordinates with two multiplications per intersection. Press Drive and turn on the voice.

The route is 3.8 km and 39 stretches of street between intersections, yet there are only ten turns: at the other intersections the road goes straight or nearly straight. The longest stretch without turns is 1837 metres. To tell a turn from a gentle bend, the program needs to know the angle and its side. Both are given by two products of vectors.

The idea

For vectors u = (u_x, u_y) and w = (w_x, w_y) the number u_xw_y − u_yw_x is positive if w is turned to the left of u (anticlockwise), negative if to the right, and zero if the vectors lie on one line. This is the cross product (also called the perp dot product). The number u_xw_x + u_yw_y, the dot product, is positive if the angle between the vectors is acute and negative if obtuse.

For three points O, A and B take the vectors OA and OB:

If a path goes O → A → B, the sign of cross(O, A, B) answers which way it turns at point A. Move the points and watch both numbers change.

The cross product also has a visual meaning: it is the signed area of the parallelogram spanned by the vectors OA and OB. When B crosses the line OA, the area passes through zero and changes sign, and that is all the “magic” there is. Half its absolute value is the area of triangle OAB, and the polygon area formula stands on this.

Why it works

Let θ be the angle through which u has to be turned anticlockwise to point the same way as w. Then

The sine is positive for a left turn through an angle from 0 to 180° and negative for a right turn. The cosine is positive while the angle is acute. The lengths of the vectors only scale the answer, so they need not be known for the sign, and no roots need be taken.

From here also comes the formula for the turn angle that the navigator uses: θ = atan2(u × w, u · w). The function atan2, given two numbers proportional to the sine and cosine, returns an angle from −180° to 180° with the right sign at once. This is better than the arccos(u · w / (|u||w|)) familiar from school: the arccosine does not know the side of the turn and loses precision near 0° and 180°, where its graph is almost vertical.

In the navigator u is the direction of entry into an intersection and w the direction of exit. If |θ| is less than 25°, the program says nothing: that is a bend in the street, not a turn. If it is more than 120°, it says “sharp”. The thresholds were chosen by eye; real navigators also take into account how many roads meet at the intersection and what they are called.

Do the segments cross

The most frequent use of orientation is checking whether segments AB and CD cross. Does the bullet cross the wall? Does the new road cross the river? Did the cursor hit a line in the drawing?

The segments cross if C and D lie on different sides of line AB, and A and B on different sides of line CD. “On different sides” means different signs of cross(A, B, C) and cross(A, B, D). One check is not enough: C and D may be on different sides of line AB, but the segment AB itself may end before CD reaches it.

What remains are the special cases on which most first solutions go wrong: one of the four products is zero. Then a point lies on the line of the other segment, and you have to check whether it lies on the segment itself, inside the rectangle spanned by its ends. Segments that touch at their ends, that lie on one line and overlap, or that lie on one line with a gap: all these cases are caught by exactly this check. Put the ends exactly on top of each other in the scene and check.

Code

The whole computational core is two one-line functions. The rest is their applications: a turn, the area of a triangle, an acute angle, a point on a segment and the intersection of segments.

The coordinates here are integers, and that is no accident. The cross product is multiplications and subtractions, with no division. For integer coordinates the answer is exact, and the sign never lies. In Go, for coordinates up to 10⁹, int64 is enough: each difference is at most 2·10⁹, a product 4·10¹⁸, the difference of two products 8·10¹⁸, and the limit of int64 is about 9.2·10¹⁸. Python does not have to think about it: its integers do not overflow.

And this is the navigator itself: Dijkstra finds a route over the street graph, the route becomes a polyline, and at each intersection atan2 of the two products decides whether to say anything to the driver.

The trap: fractions and floating point

Maps store coordinates in degrees, fractional numbers. Take three points on the line y = x: p = (0.5, 0.5), q = (12, 12), r = (24, 24). Shift p by a tiny step of 2⁻⁵³, one unit in the last place of the number 0.5 in the double format, and ask which side of line qr it ended up on. Do this for a 256 × 256 grid of shifts and compare with the exact answer.

Only 256 points lie exactly on the line, the diagonal of the grid. The computation in double is wrong for 11,972 points out of 65,536, almost one in five. Most of the errors are “on the line” where the point is slightly to one side. But for 672 points the answer is worse: the left side is called right or the other way round. The answer itself is not random, it is reproducible, but near the line it cannot be trusted.

Go has one more subtlety. The specification allows the compiler to fuse a multiplication with an addition into one FMA instruction, which rounds once instead of twice. On ARM processors the compiler does just that: the same expression without float64(…) wrappers around the products gives a different picture on ARM, 4911 errors, 4592 of them with the side mixed up. On x86 with default build settings it is the same 11,972 as in Python. The explicit conversion float64(…) forbids the fusion, and the program answers the same everywhere.

What is done about this in practice:

• Move to integers. Degrees are multiplied by 10⁷ and rounded, which gives a step of about a centimetre on the ground, finer than any GPS. This is how OpenStreetMap stores coordinates. From then on all the geometry is integer and exact.
• Compare with a tolerance. Count a point as “on the line” if |cross| is less than ε multiplied by the lengths of the vectors. This honestly says “don’t know” instead of a wrong sign, but ε has to be tuned to the problem.
• Exact predicates. First compute in double and estimate the error. If the answer is far from zero, it can be trusted; if close, recompute exactly. This is how the predicates of Jonathan Shewchuk (1997) and their descendants in geometry libraries work: almost always fast, always right.

What it costs

One orientation test is two multiplications and a few subtractions, O(1). A test for the intersection of two segments is four such tests. For a route of n segments the navigator needs O(n) operations.

It gets expensive when there are many segments. Checking all pairs of n segments is O(n²): for a hundred thousand streets that is five billion pairs. The Bentley–Ottmann sweep line algorithm finds all k intersections in O((n + k) log n), and inside it there is the same orientation test.

Where else orientation works

• Convex hull: a point is thrown off the stack if the path through it turns right.
• Point in polygon: a ray from the point crosses the sides, and how many times is counted by the segment intersection test.
• Polygon area: the sum of the cross products of neighbouring vertices, divided by two. Its sign also says whether the vertices go anticlockwise.
• Going round a polygon: a polygon is convex if all the turns along the way are to the same side.

Problems

The first problem is segment intersection with all the nasty cases: touching, overlapping, degenerate point segments and coordinates up to a billion. The second is navigator prompts on real routes, but without a single floating-point function: only the signs and comparisons of two products.