Calculate correct unit vectors for a Local Coordinate System

Question

I'm trying to implement unit vectors for a local coordinate system (LCS) based on a reference black-box implementation. In my result, I get the three unit vectors U, V, and W. The vectors are correct, but their orientation different from the referece black-box implementation.

• Input: 3-points that define the LCS. An origin, a point on the Z-axis, and a point along the XZ-plane. All given with reference to the global coordinates.

For example, for following LCS:

origin=(45.795510, -28.358400, 1.294882)
zPoint=(45.791520, -28.351800, 1.394584)
xzPoint=(45.846980, -28.443800, 1.302594)

Based on this I calculate the following U, V, W vectors:

U => (-0.039900, 0.066000, 0.997022)
V => (0.514653, -0.853923, 0.077123) 
W => (0.856469, 0.516198, 0.000104)

My calculation:

U = normalize( subtract(zPoint, origin) )
W = normalize( U x subtract(xzPoint, origin) )
V = normalize( W x U )

normalize, subtract, and x (cross-product) are same as standard vector operations.

The reference black-box implementation yields the following for the same input:

U => ( 0.514653, -0.853922, 0.077123) # V in my implementation
V => ( 0.856469,  0.516197, 0.000104) # W in my implementation
W => (-0.039900,  0.066000, 0.997021) # U in my implementation

• The problem: My implementation matches the reference in almost all cases except for the LCS above. I can't figure out what's incorrect, or if there's something missing in the input that's needed to orient these unit vectors in the right direction

Please excuse my limited math skills. Thanks a lot for all your answers

Answer 1

Mathematically, the LCS has unit vectors in the directions of the positive $X,$ $Y,$ and $Z$ axes. Let's call those unit vectors $\hat x,$ $\hat y,$ and $\hat z$ respectively.

Now zPoint is a point on the $Z$ axis of the LCS; presumably it is on the positive $Z$ axis, otherwise you cannot know which way the $Z$ axis goes. Then subtract(zPoint, origin) presumably subtracts the coordinates of origin from the coordinates of zPoint to give you the coordinates of a vector pointing from origin to zPoint, that is, in the positive $Z$ direction. Then normalize( subtract(zPoint, origin) ) makes that into a unit vector in the positive $Z$ direction, namely, $\hat z.$ So your first line of code sets U to $\hat z.$

Similarly, subtract(xzPoint, origin) gives you a vector $v$ from origin to xzPoint. Since xzPoint is a point in the $X,Z$ plane of the LCS, the vector $v$ must be $a\hat x + b\hat z$ for some real numbers $a$ and $b.$ Since you set U to $\hat z$, U x subtract(xzPoint, origin) is $$ \hat z \times v = \hat z \times (a\hat x + b\hat z) = a \hat z\times \hat x = a\hat y. $$ Then normalize( U x subtract(xzPoint, origin) ) will return $\hat y$ if $a > 0$ and $-\hat y$ if $a < 0.$ What it does if $a = 0$ is up to the function designer; it could return $0$, throw an exception, or something else. But $0\hat y = 0$ and there is no way to reliably "normalize" the zero vector to the unit $Y$ vector of your arbitrary LCS, so we had better not have $a = 0.$

So the second line of code will produce $\hat y$ in the happy circumstances where xzPoint is in the $X,Z$ plane and on the positive-$X$ side of the $Z$ axis. If we are given a point on the other side of the $Z$ axis we get a vector in the opposite direction from $\hat y$ and if xzPoint is on the $Z$ axis (which is part of the $X,Z$ plane) then we don't get an axis vector at all. But whatever you get from the second line, you set W to it.

Assuming we were given a "happy" choice of xzPoint, we now have U set to $\hat z$ and W set to $\hat y.$ Then the third line computes $$\hat y \times \hat z = \hat x.$$ Note that in the absence of numerical error, this equation is exactly correct; the cross product of two orthogonal unit vectors is a unit vector. But you might have some numerical error in the calculations of W (aka $\hat y$) and U (aka $\hat z$), so you normalize the result of the cross product in order to have a more accurate value of $\hat x.$ You then set V to $\hat x$ in the third line of code.

So now you have:

• V set to $\hat x.$
• W set to $\hat y.$
• U set to $\hat z.$

The reference implementation, on the other hand, sets U to $\hat x,$ V to $\hat y,$ and W to $\hat z,$ which we might expect since this preserves the alphabetical order of the vectors relative to the axes.

The main discrepancy is how you named your vectors U, V, W. If you want to match the reference implementation, you want W to be $\hat z,$ so your first line of code should set W, not U, and so forth. I also see some differences in the sixth decimal place, possibly due to different algorithms with different numerical errors.