Freesteel Blog » 2006 » November

Wednesday, November 1st, 2006 at 5:46 pm - - Machining

I got sent a bug which took two days to track down through a set of too-optimized code that I’d hoped was so solid I’d never need to look at it again — until I have to implement new cutter shapes, like tapered cutters. (Maybe now is a good time to do it. Why cares; I’m not in the mood.)

The problem has to do with testing the tip of the cutter against a single triangle. In this case, the triangle has one edge that has a Z-displacement of 3e-15, due to something less than perfect happening in the triangulator. This is an example of a not-quite-perfectly-horizontal edge of the triangle. Marginal values like these are very good at breaking calculations. Ordinary edges don’t cause problems. Perfectly aligned edges get treated as a special case. It’s the ones that aren’t quite perfectly aligned that are difficult when the algorithm you put it through gets into trouble with the precision.

In geometric algorithms, there are always a quick and dirty methods for dealing with problems using epsilon values; you pick a small number, like 1e-8, and say that when two Z-values are closer than this, the line is horizontal, and treate it accordingly.

A similar method is to filter all the points that come in by rounding all coordinates to the 5th decimal place. The boys in Denmark put that in, but I keep it turned off in my versions. It’s beyond the precision of any machine tool, but it’s the principle that counts.

I don’t like these epsilon tricks, even though they result in apparently robust code with the minimum of effort. It’s a bad habit. The error is showing a weakness in the algorithm. Usually the filter gets put in by somebody who doesn’t have the time or interest to investigate the weakness in some obscure and complex algorithm, and what they’re doing is changing the inputs to displace it away from the weakness. That’s not good in the long term, because the fix gets used by everything else. The next geometric algorithm to be worked on will see this filtered data and can get away with being more sloppy because the input data doesn’t now test it properly. Maybe it’s so shoddy that that the flaws in it show through until they get papered over by revising epsilon upwards slightly to get it away from the difficult areas. And so on, until the rot really sets in. At best you have calculations that run differently in inches or millimetres. At worst, you end up patching over problems on a case by case basis for decades unable to recognize that it’s all due to a tiny flaw in the calculation at the very start of the development that should have been fixed 20 years ago.

It’s better to refrain from using epsilon values in any geometric calculations. They’re often not needed, though they tend to be the tool of first choice of programmers who don’t own their code and don’t really care. It’s something to look for on inspection, if you happen to suddenly gain access to a large body of machining code, like if you come in to a company as a new employee.

The triangle in question

I have points a,b,c of a triangle. It’s critical that I know the orientation of the triangle projected to the planes perpendicular to all three axes, for reasons of knowing which way to offset the edges to make them go outwards. The Z-orientation refers to whether the path abca goes clockwise or counter-clockwise, as seen in the XY plane. If it’s counter-clockwise, then, taking the line from a to b, turning right (in the XY plane) takes me to the outside of the triangle. I need these orientations for the YZ-plane projection and the ZX-plane projection (note that the XZ-plane is a mirroring, because XZ is not in the right order within the cycle XYZX).

The solution is elementary. Let n=(b - a)x(c - a), the cross-product of the two vectors along two sides of the triangle. The vector n is normal to the triangle, and, if a,b,c all lie in a line and do not form a triangle, it has length zero.

If n_z is positive, then the normal points upwards, and abca goes counter-clockwise in the XY plane. If n_z is negative, it’s clockwise. If n_z is zero, then the three points are colinear in the XY plane (it’s aligned vertically).

The same is true for the other two axis. This is an effective way of working out their orientations.

The problem I got was when n_x became zero for a particular triangle, and the orientation was wrong, because I treated the zero value as though it was positive.

If we look at the equation for the X-coordinate of the normal vector, it’s

    n_x = (b_y - a_y)*(c_z - a_z) - (b_z - a_z)*(c_y - a_y) 

This expands to

    n_x = b_y*c_z - b_y*a_z - a_y*c_z - b_z*c_y + b_z*a_y + a_z*c_y 

because the a_y*a_z cancels out. This is a good thing, because it removes one source of imprecision.

Suppose the values of the vectors in the YZ plane are:

    a_yz=(0, 10), b_yz=(10, 0), c_yz=(10, 1e-15)

Clearly, this is a counter-clockwise triangle, because c is higher than b along the same constant Y-line, and a is to the left.

Substituting into the equation gives:

    n_x = 1e-14 - 100 - 0 - 0 + 0 + 100

If you are unlucky, the (1e-14 - 100) will give -100 exactly, due to a rounding error, and adding 100 to that gives 0, which doesn’t tell you which way the triangle is oriented in the YZ plane, and there’s a 50% chance you’ll get it wrong, and the cutter won’t know it should have hit the triangle, and the resulting cutter-path will gouge the model.

I’ve added an extra function to work out the sum properly, making

    n_x = Sum6Acc(b_y*c_z, -b_y*a_z, -a_y*c_z, -b_z*c_y, b_z*a_y, a_z*c_y) 

where:

Sum6Acc(a[0], a[1], a[2], a[3], a[4], a[5])
{
    sort(a[], abs())  // by ascending absolute value
    return (((((a[5] + a[4]) + a[3]) + a[2]) + a[1]) + a[0]) 
}

This allows the large values to cancel out before they swamp the smaller terms. Luckily, the triangle normals are calculated once and stored in memory, so we don’t have to run through this for each triangle each time a cutter shape is tested against it.

It’s a result. I’m not completely happy with it because it still looks a potentially weak. At least there are no epsilons, and I’ve now got a few more ASSERTS in place to detect instances of misorientation when they come through. Before now there wasn’t enough, which was why the problem didn’t show immediately when I ran it in debug mode. That was the really annoying part of it. I found at least one ASSERT that was pointing it out, before I weakened it because it was firing too much.