Freesteel Blog » Triangle orientation interlude

Triangle orientation interlude

Saturday, May 3rd, 2008 at 1:12 am Written by:

It’s critical to be able to calculate the orientation of a triangle correctly, particularly if it’s the long thin kind where the endpoints differ by 1e-15 in the y direction. The orientation is either clockwise or anticlockwise, and it should be the same no matter which of the three vertices you start with.

For a triangle with points a, b and c, we consider the vectors u = b - a and v = c - a. If

CPerp(u) . v

is positive, then the triangle is oriented clockwise, where CPerp((x,y)) = (y, -x) is the function that turns a 2-dimensional vector 90 degrees clockwise.

Writing this in individual coordinates gives:

CPerp(b - a) . (c - a)
    = (by - ay, ax - bx) . (cx - ax, cy - ay)
    = (by - ay) (cx - ax) + (ax - bx) (cy - ay)

which you can either expand out into 8 terms (two of which are ±ayax which cancel immediately) and add them up, or you can evaluate them as they stand. You can get different answers; sometimes one is zero and the other is non-zero, and since the sign of the result determins the orientation of the triangle, we pick the non-zero answer.

To reduce loss of precision, sumations are done from large to small; after putting the six terms into the array double s[6];I write:

sort(s, s + 6, FabsOrder());
double sum = (((((s[5] + s[4]) + s[3]) + s[2]) + s[1]) + s[0]);

where

struct FabsOrder {
bool operator()(double a, double b)
{ return fabs(a) < fabs(b); } };

That’s the C++ done with. Now the problem. There was a conflict in the sign between one value and the other which was highlighted by the debug code, but, because it seemed infinitesimally small and we’d forgotten the significance, we commented it out and encountered a situation much later in the algorithm where there was a gouge that took a while to isolate.

The calculation problem.

In the above equation, if by = cy, then the the alternative ways to calculate the result are given by the two sides of the equation:

p (r - s) - q (r - s) = p r - p s - q r + q s

And if r - s is positive, then the result should be positive whenever p > q

But here’s a table of calculations that happens.

s 13.014700000000001 * t – 13.014699999999999 * t
20.0 5.6843418860808015e-14
20.677 0.0
21.0 5.6843418860808015e-14
1e15 2.0

Now I had been counting on the right hand column to be increasing, because 13.014700000000001 - 13.014699999999999 > 0, but unfortunately it does not. This explains how I can get the wrong sign in the answer, rather than 0.0 when the precision falls off the end, because I am taking the difference of two values where s is increasing in the expectation that the difference of the products will increase too.

This has been a day’s work. I’ve marked this down as yet another of many things to watch out for with floating point arithmetic.

Meanwhile, after scratching my head all day over this problem, the BBC reported that some computer keyboards are dirtier than a toilet seat. If I saw a toilet seat as grubby as this keyboard is now, I wouldn’t sit on it.

Leave a comment

XHTML: You can use these tags: <a href="" title=""> <abbr title=""> <acronym title=""> <blockquote cite=""> <code> <em> <strong>