Month: June 2014

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.

Whether p is above, exactly on, or below line p, r can be determined from the orientation of triangle p, q, r.

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 floats. 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 floats 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.