Rob Smallshire from Good With Computers

Computational geometry - a world where lines have zero thickness,
circles are perfectly round and points are dimensionless. Creating
robust geometric algorithms using finite precision number types such as
`float` is fiendishly difficult because it's not possible to exactly
represent numbers such as one-third, which rather gets in the way of
performing seemingly simple operations like dividing a line into exactly
three equal segments. In this short series of posts, we'll look at some
of the pitfalls of geometric computation, with examples in Python,
although the key messages are true with finite-precision floating point
numbers in any language.

Rational numbers, modelled by Python's `Fraction` [1] type can be useful for
implementing *robust* geometric algorithms. These algorithms are often deeply
elegant and surprising because they must avoid any detour into the realm of
irrational numbers which cannot be represented in finite precision, which means
that using seemingly innocuous operations like square root, for example to
determine the length of a line using Pythagoras, are not permitted.

One algorithm which benefits from rational numbers is a simple collinearity test
determining whether three points lie on the same line. This can be further
refined to consider whether a query point \(p\) is above, exactly on,
or below the line. Now there are many ways to approach writing such a function,
and like many questions in computational geometry the naïve approaches are
either overly complex, inefficient, or just plain wrong, albeit in subtle ways.
I won't cover the story of how to arrive at a robust algorithm, that story is
entertaining covered in *Writing Programs for "The Book"* by Brian Hayes. [2]
Rather, we'll start where Brian leaves off by showing how to implement the
algorithm in Python using both floating-point and exact arithmetic so we can
understand the tradeoffs between performance and correctness inherent in these
choices. Along the way, we'll perhaps touch on some aspects of Python which may
be new to you.

You don't need to understand the mathematics of the orientation test to appreciate the point of what we're about to demonstrate, suffice to say that the orientation of three two-dimensional points can be concisely found by computing the sign of the determinant of a three by three matrix containing the \(x\) and \(y\) coordinates of the points in question, where the determinant happens to be the signed area of the triangle formed by the three points:

The function \(o\) returns \(+1\) if the polyline \(p\), \(q\), \(r\) executes a left turn and the loop is counterclockwise, \(0\) if the polyline is straight, or \(-1\) if the polyline executes a right turn and the loop is clockwise. These values can in turn be interpreted in terms of whether the query point \(p\) is above, on, or below the line through \(q\) and \(r\).

To cast this formula in Python, we need a sign function and a means of
computing the determinant. Both of these are straightforward, although
perhaps not obvious, and give us the opportunity to explore a little
appreciated aspect of Python. First, the `sign()` function. You may be
surprised to learn − and you wouldn't be alone − that there is no
built-in or library function in Python which returns the sign of a
number as \(-1\), \(0\) or \(+1\), so
we need to roll our own. The simplest solution is probably something
like this:

```
>>> def sign(x):
... if x < 0:
... return -1
... elif x > 0:
... return 1
... return 0
...
>>> sign(5)
1
>>> sign(-5)
-1
>>> sign(0)
0
```

This works well enough. A more elegant solution would be to exploit an
interesting behaviour of the `bool` type, specifically how it behaves
under subtraction. Let's do a few experiments:

```
>>> False - False
0
>>> False - True
-1
>>> True - False
1
>>> True - True
0
```

Intriguingly, subtraction of `bool` objects has an integer result! In
fact, when used in arithmetic operations this way, `True` is
equivalent to positive one and `False` is equivalent to zero. We can
use this behaviour to implement a most elegant `sign()` function:

```
>>> def sign(x):
... return (x > 0) - (x < 0)
...
>>> sign(-5)
-1
>>> sign(5)
1
>>> sign(0)
0
```

Now we need to compute the determinant. In our case this turns out to reduce down to simply:

so the definition of our `orientation()` function using tuple
coordinate pairs for each point becomes just:

```
def orientation(p, q, r):
d = (q[0] - p[0]) * (r[1] - p[1]) - (q[1] - p[1]) * (r[0] - p[0])
return sign(d)
```

Let's test this on on some examples. First we set up three points `a`,
`b` and `c`:

```
>>> a = (0, 0)
>>> b = (4, 0)
>>> c = (4, 3)
```

Now we test the orientation of a ➔ b ➔ c:

```
>>> orientation(a, b, c)
1
```

This represents a left turn, so the function returns positive one. On the other hand the orientation of a ➔ c ➔ b is negative one:

```
>>> orientation(a, c, b)
-1
```

Let's introduce a fourth point `d` which is collinear with `a` and
`c`. As expected our `orientation()` function returns zero for the
group a ➔ c ➔ d:

```
>>> d = (8, 6)
>>> orientation(a, c, d)
0
```

So far so good. Everything we have done so far is using integer numbers
which, in Python, have arbitrary precision. Our function only uses
multiplication and subtraction, with no division to result in `float`
values, so all of that precision is preserved. But what happens if we
use floating point values as our input data? Let's try some different
values using `float`s. Here are three points which lie on a diagonal
line:

```
>>> e = (0.5, 0.5)
>>> f = (12.0, 12.0)
>>> g = (24.0, 24.0)
```

As we would expect, our orientation test determines that these points are collinear:

```
>>> orientation(e, f, g)
0
```

Furthermore, moving the point `e` up a little, by increasing its
\(y\) coordinate by even a tiny amount, gives the answer we
would expect:

```
>>> e = (0.5, 0.5000000000000018)
>>> orientation(e, f, g)
1
```

Now let's increase the \(y\) coordinate just a little more. In fact, we'll increase it by the smallest possible amount to the next representable [3] floating point number:

```
>>> e = (0.5, 0.5000000000000019)
>>> orientation(e, f, g)
0
```

Wow! According to our orientation function the points `e`, `f` and
`g` are collinear again. This cannot possibly be! In fact, we can go
through the next 23 successive floating point values up to,

```
>>> e = (0.5, 0.5000000000000044)
>>> orientation(e, f, g)
0
```

with our function still reporting that the three points are collinear, until we get to this value,

```
>>> e = (0.5, 0.5000000000000046)
```

at which point things settle down and become well behaved again:

```
>>> orientation(e, f, g)
1
```

What's happening here is that we've run into problems with the finite
precision of Python `float`s at points very close the diagonal line,
and the mathematical assumptions we make in our formula about how
numbers work break down due to the fact that floating point numbers are
a far from a perfect model of real numbers. Next time, we'll investigate
more thoroughly, the behaviour of this orientation test at the limits of
floating-point precision.

[1] | The Fraction type which models rational numbers is defined in the Python Standard Library fractions module. |

[2] | Hayes, Brian. (2007) Writing Programs for "The Book". In: Oram, A. & Wilson, G., eds. Beautiful Code O'Reilly Media. pp. 539–551. |

[3] | In C we could use the nextafter() function to generate the next
representable floating point number. Unfortunately, nextafter() is
not exposed to Python. Various workarounds are available, including a version built
into Numpy, directly calling the C version using ctypes and a pure
Python implementation. |