Safe vector angle calculations

Discuss programming topics for any language, any source base. If it is programming related but doesn't fit in one of the below categories, it goes here.
Post Reply
taniwha
Posts: 400
Joined: Thu Jan 14, 2010 7:11 am
Contact:

Safe vector angle calculations

Post by taniwha »

I'm sure we've all been bitten by acos(x) or asin(x) producing nan (or an exception) despite efforts to ensure x is in the -1..1 range. I learned a very nice bit of math while I was off in KSP land. Don't use acos or asin. Use atan, or even better, atan2. (Indirect thanks to William Kahan, more direct thanks to "egg" (KSP modder (Principia mod)).

For two unit vectors a and b:

Code: Select all

ang = 2 * atan2 (mag(a - b), mag(a + b);// use 360/pi for degrees
(mag(v) is the length of v)

This is perfectly safe at any angle, and can never produce a nan unless either a or b has a nan in a component.

If either a or b may be a non-unit vector (even zero!!!), then with a slight modification the angle can still be safely determined:

Code: Select all

ang = 2 * atan2(mag(mag(b)*a - mag(a)*b), mag(mag(b)*a + mag(a)*b));
Of course, the obvious optimization of caching mag(b)*a and mag(a)*b can be quite helpful.

There are only two disadvantages to using this: it's a little slower (atan and thus atan2 (which is just a very smart wrapper for atan) is inherently slower than acos or asin), and you lose the direction of the angle. However, this can be retrieved by finding the cross-product of a and b and dotting that with a reference direction to get the sign of the angle.

Some nice bonus information...
For unit vectors, mag(a - b) == 2*sin(ang(a,b)/2), ie, twice the sin of the half-angle between a and b, and mag(a + b) == 2*cos(ang(a,b)/2) (twice the cos of the half-angle)

Now, this is very handy when it comes to finding the quaternion describing the rotation between the two vectors. When I implemented the TINU (There Is No Up) mod for KSP, I discovered that Unity's Quaternion.FromToRotation was not usable for very small angles (think per-frame rotation while in low orbit, worse in high orbit). The reason for this is, I assume, due to using epsilons and guarding against going outside the -1..1 range for acos.

I'll leave the derivation as an exercise (hint, it takes advantage of trig identities), but this is taken directly from the code in TINU:

Code: Select all

Quaternion fromtorot(Vector3 a, Vector3 b)
{
	float ma = a.magnitude;
	float mb = b.magnitude;
	Vector3 mb_a = mb * a;
	Vector3 ma_b = ma * b;
	float den = 2 * ma * mb;
	float mba_mab = (mb_a + ma_b).magnitude;
	float c = mba_mab / den;
	Vector3 v = Vector3.Cross (a, b) / mba_mab;
	return new Quaternion(v.x, v.y, v.z, c);
}
Only an exact 180 degree rotation (ie, a and b are anti-parallel) or either vector is invalid (nan, inf, 0) can produce a bogus quaternion. Also, the code is chirality agnostic.
Leave others their otherness.
http://quakeforge.net/
Post Reply