9

What is the best (fastest) way to compute two vectors that are perpendicular to the third vector(X) and also perpendicular to each other?

This is how am I computing this vectors right now:

// HELPER - unit vector that is NOT parallel to X
x_axis = normalize(X);
y_axis = crossProduct(x_axis, HELPER);
z_axis = crossProduct(x_axis, y_axis);

I know there is infinite number of solutions to this, and I don't care which one will be my solution.

What is behind this question: I need to construct transformation matrix, where I know which direction should X axis (first column in matrix) be pointing. I need to calculate Y and Z axis (second and third column). As we know, all axes must be perpendicular to each other.

sidon
  • 1,434
  • 1
  • 17
  • 30
  • Do any of the algorithms [here](http://blog.selfshadow.com/2011/10/17/perp-vectors/) have what you need? They concern generating the y_axis in a robust and efficient way. – John McFarlane Mar 26 '13 at 01:08
  • This is a almost the same question as: [http://stackoverflow.com/questions/19649452/given-a-single-arbitrary-unit-vector-what-is-the-best-method-to-compute-an-arbi/22438742#22438742][this] with one ortho vector, just cross product to produce the second. – ideasman42 Mar 16 '14 at 15:27

5 Answers5

7

What I have done, provided that X<>0 or Y<>0 is

  1. A = [-Y, X, 0]
  2. B = [-X*Z, -Y*Z, X*X+Y*Y]

and then normalize the vectors.

[ X,Y,Z]·[-Y,X,0] = -X*Y+Y*X = 0
[ X,Y,Z]·[-X*Z,-Y*Z,X*X+Y*Y] = -X*X*Z-Y*Y*Z+Z*(X*X+Y*Y) = 0
[-Y,X,0]·[-X*Z,-Y*Z,X*X+Y*Y] = Y*X*Z+X*Y*Z = 0

This is called the nullspace of your vector.

If X=0 and Y=0 then A=[1,0,0], B=[0,1,0].

John Alexiou
  • 28,472
  • 11
  • 77
  • 133
  • And yet another mathematically exact, but numerically unstable computation gets introduced. –  May 22 '12 at 18:20
  • Please explain "unstable". Can you provide with an (non-zero) example were it fails? – John Alexiou May 22 '12 at 19:53
  • For example, if X and Y are different in their magnitudes, squaring them will tend to yield results that while algebraically correct, are not numerically correct. X*X + Y*Y reduces to Y*Y when Y is reasonably large enough compared to X. –  May 22 '12 at 21:55
  • The input vector is normalized according to the OP, and so the worst case would be `[1, ε, 0]` where ε is the machine precision. The results would be `[-ε,1,0]` and `[0,0,1]` which is still the _nullspace_. I really don't see your point here at all (in this case). – John Alexiou May 23 '12 at 00:47
  • You might think the the `[ε,ε,1]` case fails, but when the results are normalized in the last step it works out correctly. Don;t forger modern computers compute with 80-bits floats internally, but only display 32 or 64 bits in the end. – John Alexiou May 23 '12 at 00:55
3

This is the way to do it.
It's also probably the only way to do it. Any other way would be mathematically equivalent.
It may be possible to save a few cycles by opening the crossProduct computation and making sure you're not doing the same multiplications more than once but that's really far into micro-optimization land.

One thing you should be careful is of course the HELPER vector. Not only does it has to be not parallel to X but it's also a good idea that it would be VERY not parallel to X. If X and HELPER are going to be even somewhat parallel, your floating point calculation is going to be unstable and inaccurate. You can test and see what happens if the dot product of X and HELPER is something like 0.9999.

shoosh
  • 76,898
  • 55
  • 205
  • 325
  • Thanks for your response. Yes HELPER is nowhere near to be parallel to X. I'm developing for android platform and I was hoping that there would be some trick to compute at least the first parallel vector somewhat cheaper. I'll give this question some time and if no better answer will be made, i'll accept this as an answer. – sidon May 22 '12 at 12:56
  • No. It is NOT the ONLY way to do it, nor is it necessarily the best way, depending on the value of HELPER. –  May 22 '12 at 18:16
3

There is a method to find a good HELPER (really - it is ready to be your y_axis).

Let's X = (ax, ay, az). Choose 2 elements with bigger magnitude, exchange them, and negate one of them. Set to zero third element (with the least magnitude). This vector is perpendicular to X.

Example:

if (ax <= ay) and (ax <= az) then HELPER = (0, -az, ay) (or (0, az, -ay))

X*HELPER = 0*0 - ay*az + az*ay = 0

if (ay <= ax) and (ay <= az) then HELPER = (az, 0, -ay)

MBo
  • 77,366
  • 5
  • 53
  • 86
1

For a good HELPER vector: find the coordinate of X with the smallest absolute value, and use that coordinate axis:

absX = abs(X.x); absY = abs(X.y); absZ = abs(X.z);
if(absX < absY) {
  if(absZ < absX)
    HELPER = vector(0,0,1);
  else // absX <= absZ
    HELPER = vector(1,0,0);
} else { // absY <= absX
  if(absZ < absY)
    HELPER = vector(0,0,1);
  else // absY <= absZ
    HELPER = vector(0,1,0);
}

Note: this is effectively very similar to @MBo's answer: taking the cross-product with the smallest coordinate axis is equivalent to setting the smallest coordinate to zero, exchanging the larger two, and negating one.

comingstorm
  • 25,557
  • 3
  • 43
  • 67
0

I think the minimum maximum magnatude out of all element in a unit vector is always greater than 0.577, so you may be able to get away with this:

-> Reduce the problem of finding a perpendicular vector to a 3D vector to a 2D vector by finding any element whose magnatude is greater than say 0.5, then ignore a different element (use 0 in its place) and apply the perpendicular to a 2D vector formula in the remaining elements (for 2D x-axis=(ax,ay) -> y-axis=(-ay,ax))

let x-axis be represented by (ax,ay,az)

if (abs(ay) > 0.5) {
  y-axis = normalize((-ay,ax,0))
} else if (abs(az) > 0.5) {
  y-axis = normalize((0,-az,ay))
} else if (abs(ax) > 0.5) {
  y-axis = normalize((az,0,-ax))
} else {
  error("Impossible unit vector")
}
clinux
  • 2,984
  • 2
  • 23
  • 26