Freesteel Blog » Perpendicular unit vector frame in 3D

Perpendicular unit vector frame in 3D

Tuesday, March 26th, 2013 at 9:44 pm Written by:

Let’s begin with what a unit vector is, and let’s not get sloppy. If we can’t be bothered to do this right, then what does it say about our other more advanced geometric algorithms?

I’m a computational geometry programmer, and I don’t like to type, so my 3D point class is defined like so:

class P3 { double x, y, z }; 

I don’t think much of a special type called vector3 that is exactly the same thing, except you can add two of vector3s together, but you can’t add two P3s together. This sort of excessive type separation is annoying, serves no useful purpose, and is often promoted by the same folk who think that any attention to the detail of making correct unit vectors is not important — even though it actually does make a difference to the end calculation as well as the functions you might use to separate and optimize different cases.

If n is a P3 point and n is a unit vector, then abs(Len(n) – 1) < epsilon (ie it is approximately equal to 1).

We can also assume that 1 >= abs(n.x), abs(n.y), abs(n.z) should be true.

Here is the contentious part: if n.x = n.y = 0 then n.z = +-1.

It’s not 1.000000001 or 0.999999998, it’s exactly 1 when you have a cardinal vector. It’s a special case, so you can then use the condition (n == P3(0,0,1)) to separate it out and optimize accordingly by avoiding doing the transform at all because your vector is already aligned.

How do these shoddy unit vectors get into the system in the first place? Usually through the shoddy normalizing function:

P3 Normalize(P3 p)  { return p / Len(p); }

Actually, it’s the Len(p) which is the problem, because this is normally implemented as:

double Len(P3 p) { return sqrt(p.x*p.x + p.y*p.y + p.z*p.z); }

so if p = P3(0, 0, 3.432453452352348844), then p.z is squared and then square rooted, which means it can end up on either side of the original value, which, when divided into the original, means you can get something slightly greater or slightly less than 1.

It would be more useful for the Len(p) function to recognize the cardinal vectors (when it matters, as does in the Normalize function), and return exactly abs(p.z) when p.x = p.y = 0. There is no reason not to assume that floating point arithmetic will give us the value 1 when a number is divided by itself.

But this wasn’t even where the problem was. Attention to this detail revealed a bug in my PerpVecsN(tz) function that I use to create a coordinate space from the position of a 5-axis tool.

The tool is rotationally symmetric (because it is a milling piece), and I am given the vector tz, which is the direction of the axis of the tool, but all my milling tool collision functions are built for 3-axis where the tool is assumed to be parallel to the z axis.

So, what I need are two vectors, tx and ty that are orthogonal to each other and to tz to use as a coordinate space. Then the transform I will apply to a point of the model p to move it to the position where my 3-axis collision functions will work is P3(Dot(p, tx), Dot(p, ty), Dot(p, tz)).

(This is the closest I get to allowing matrices in my geometric code. Once you allow those beasts in as general objects, you don’t know what they’re doing. Here I know it’s an orthogonal transform, and I understand it well enough to develop culling routines on regions of the model that I can avoid transforming in the first place when the collision tool is nowhere near it. That’s the sort of thing to keep in mind with all these functions — the fastest calculations are the ones you work out a way to not do.)

The function I made (assuming that tz is a unit vector) was:

pair<P3, P3> PerpVecsN(P3 tz)
{
    P3 tx = Normalize(abs(tz.z) > 0.3 ? P3(tz.z, 0.0, -tz.x) : P3(tz.y, -tz.x, 0.0));
    P3 ty = CrossProd(tz, tx);
    return pair<P3, P3>(tx, ty);
}

This works by quickly finding a perpendicular vector in the XY or XZ planes, depending on the alignment of the incoming tz, and then calculating the perpendicular vector to them by knowing that the cross product of two unit length perpendicular vectors is itself a unit vector.

This means I get to the result by a single square root and a few multiplications and additions.

But I was getting all kinds of failures with the unitization of ty, with coordinates that were like (0, 0.999999999999, 0), which I don’t like.

So I expanded the cross product to work out what was going on.

Firstly, let’s assume that abs(tz.z) > 0.3. Then:

double txl = Len(P2(tz.z, tz.x));
P3 tx = P3(tz.z / txl, 0.0, -tz.x / txl);

Obviously:

Dot(tx, tz) = (tz.z * tz.x / txl - tz.x * tz.z / txl) = 0

We can also expand the cross product between tx and tz, taking account of the fact that tx.y = 0:

P3 ty = P3(-(tz.y * tz.x) / txl, txl, (-tz.y * tz.z) / txl);

Obviously:

Dot(tx, tz) = -tz.z * tz.y * tz.x / txl^2 + tz.x * tz.y * tz.z / txl^2 = 0

Dot(ty, tz) = -(tz.y * tz.x * tz.x) / txl + txl * tz.y - (tz.y * tz.z * tz.z) / txl
            = (-tz.x^2 - tz.z^2 + txl^2) * tz.y / txl
            = (-txl^2 + txl^2) * tz.y / txl
            = 0

And the bit I was getting hung up on was the case when tz.y = 0, and my txl = Len(P2(tz.z, tz.x)) was returning 0.9999999999999 because P2(tz.z, tz.x) was a unit vector, and ty was then (0, 0.99999999999999, 0), which I don’t like.

Finally, with some help, we made tz.y = 0 into a special case. But I was losing it for a while, as the miscalculation only occurred on a different machine in completely optimized code, and it was impossible for me to recreate it.

Well, anyway, it’s a way of life. And it seems to work in that I am not continually debugging the same code for the rest of my life — though this might just be because I am lazy. I have been hunting around the Company for a while looking for someone else, anyone, with the same attitude towards computational geometry arithmetic as I have, to no avail. Surely there’s got to be some variety of opinion out there. You can’t all use epsilon fudge values without questioning it. Remember when the goto was an important part of all programming languages? This is like that. It’s a radical suggestion to make the case that it is harmful. One day, if there is any progressive development in the field (ie people learning from each other’s long term mistakes), this point of view will get a proper airing.

Leave a comment

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