Rational Orthogonal Unit Vectors
The problem of creating a correctly oriented, orthonormal basis from a single unit vector is surprisingly common in geometry/graphics code. Specifically, we have some unit vector \( v \) and wish to generate two other unit vectors \(a\) and \(b\) so that \(a\cdot b = 0\). and \(a\times b = v\). This usually happens when we have radial symmetry about an axis in 3D and don't have a canonical "orthogonal" direction. It's common in ray tracing with radially symmetric BSDFs, for example. Or, in my case, generating a 4 x 4 affine transform matrix from an axis of rotation and a rotation amount.
I've come across this problem innumerable times before, and my typical solution is to come up with any vector \(w\) that isn't colinear with \( v \), then set \(a=\frac{v\times w}{\|v\times w\|}\) and \(b=v\times a\). Here \(w\) would be one of the three basis vectors, chosen by the coordinate of \(v\) with the least magnitude.
I've always found it rather unfortunate to resort to casework for such a simple problem. And the square root required for \(\|v\times w\|\) is expensive. Of course, I'm not exactly a pioneer here; Jeppe Revall Frisvad wrote about this in 2012 in Building an Orthonormal Basis from a 3D Unit Vector Without Normalization. In the paper, Frisvad presents a neat formula for finding \(a\) and \(b\) without a square root. Reproduced in terms of my notation: \[ a=\left( 1-\frac{v_x^2}{1+v_z}, -\frac{v_xv_y}{1+v_z}, -v_x \right) \] \[ b=\left( -\frac{v_xv_y}{1+v_z}, 1-\frac{v_y^2}{1+v_z}, -v_y \right) \] He uses a quaternion formula from "Real Time Rendering" (involving a square root), does some algebra, and notes that "the square root disappears" as he reaches his slick result.
Intellectually, I find this unsatisfying, if impressive. Let's instead spend a little more time classifying all rational expressions that generate orthonormal vectors and see where Frisvad's formula lands there.
Setup
Technically speaking, we'd like to classify unit orthogonal vectors \(a\) in \(\mathbb{Q}(v_x,v_y,v_z)^3\), that is, three-vectors in the field extension of \(\mathbb{Q}\) generated by \(v_x\), \(v_y\), and \(v_z\). However, rational points on the unit sphere are dense over \(S^2\), and our floating-point representation for \(v\) approximate an infinite supply of rational unit vectors just as well as the irrational ones. Therefore, we assume without loss of generality that \(v\) has rational coordinates. Now we can ditch the field extensions and classify (unit orthogonal) vectors \(a\) in \(\mathbb{Q}^3\).
Stereographic Projection
This whole exercise is really about the stereographic projection from \(S^2\) (minus a pole) to the plane. Why? Well, the following map is a bijection(!) from 3D rational points on \(S^2\setminus \{(0,0,1)\}\) to 2D rational points on the plane: \[ (x,y,z) \mapsto \left( \frac{x}{1-z}, \frac{y}{1-z} \right) \] with inverse \[ (x,y)\mapsto \frac{1}{1+x^2+y^2}\left(2x,2y,x^2+y^2-1\right)\,. \] The fact that the projection and its inverse preserves rational coordinates means that rational points on the sphere (minus the north pole) map exactly to rational points on the plane. (Incidentally, this is an easy way to see that rational points are dense on \(S^2\).)
Lets consider rational \(v=(x,y,z)\) (not the north pole) mapping to some \((X,Y)\) on the plane via stereographic projection. A possible rational unit vector \(a=(s,t,u)\) (not on the north pole) mapping to \((S,T)\) on the plane is orthogonal to \(v\) if and only if \[ \begin{align*} v\cdot a=0 &\iff xs+yt+zu = 0 \\ &\iff \frac{1}{1+X^2+Y^2}\cdot\frac{1}{1+S^2+T^2}(4XS+4YT+(X^2+Y^2-1)(S^2+T^2-1)) = 0 \\ &\iff 4XS+4YT+(X^2+Y^2-1)(S^2+T^2-1) = 0\\ &\iff ((X^2+Y^2-1)S+2X)^2 + ((X^2+Y^2-1)T+2Y)^2 = 4X^2+4Y^2+(X^2+Y^2-1)^2 \\ &\iff ((X^2+Y^2-1)S+2X)^2 + ((X^2+Y^2-1)T+2Y)^2 = (X^2+Y^2+1)^2 \\ &\iff \left(\frac{X^2+Y^2-1}{X^2+Y^2+1}S+\frac{2X}{X^2+Y^2+1}\right)^2+\left(\frac{X^2+Y^2-1}{X^2+Y^2+1}T+\frac{2Y}{X^2+Y^2-1}\right)^2=1 \\ &\iff (zS+x)^2+(zT+y)^2=1\,. \end{align*} \] So the rational point \((zS+x,zT+y)\) lies on the unit circle. Parameterizing by an angle \(\theta\), we have \[ \begin{align*} zS+x=\cos\theta\implies S=\frac{\cos(\theta)-x}{z} \\ zT+y=\sin\theta\implies T=\frac{\sin(\theta)-y}{z} \end{align*}\,. \] Sending this back to the unit sphere, we see that \[ (s,t,u)=\frac{(2z(\cos(\theta)-x), 2z(\sin(\theta)-y), 1+x^2+y^2-z^2-2x\cos(\theta)-2y\sin(\theta))}{1+x^2+y^2+z^2-2x\cos(\theta)-2y\sin(\theta)}\,. \] Since \(v=(x,y,z)\) has length one, we have \(x^2+y^2+z^2=1\) and \[ (s,t,u)=\frac{(z(\cos(\theta)-x), z(\sin(\theta)-y), 1-z^2-x\cos\theta-y\sin\theta)}{1-x\cos\theta-y\sin\theta}\,. \]
Recall that \(\cos\theta\) and \(\sin\theta\) are rational. Pythagorean triples will generate such rational points on the unit circle, which yields a family of (all!) rational expressions for orthogonal unit vectors to \(v\). Here are a few examples:
| Pythagorean Triple | \((\cos\theta,\sin\theta)\) | \((s,t,u)\) |
|---|---|---|
| \((1,0,1)\) | \((1,0)\) | \(\left(z,-\frac{yz}{x-1},\frac{x-1+z^2}{x-1}\right)\) |
| \((0,1,1)\) | \((0,1)\) | \(\left(-\frac{xz}{y-1},z,\frac{y-1+z^2}{y-1}\right)\) |
| \((3,4,5)\) | \((3/5,4/5)\) | \(\left(\frac{z(5x-3)}{3x+4y-5}, \frac{z(5y-4)}{3x+4y-5}, \frac{3x+4y-5+5z^2}{3x+4y-5}\right)\) |
| \((5,12,13)\) | \((5/13,12/13)\) | \( \left(\frac{z(13x-5)}{5x+12y-13}, \frac{z(13y-12)}{5x+12y-13}, \frac{5x+12y-13+13z^2}{5x+12y-13}\right) \) |
The other basis vector \(b\) is, of course, only a cross product away.
Degeneracies
Observe that even with \((x,y,z)=(0,0,1)\), we get \[ (s,t,u)=(\cos\theta, \sin\theta, 0) \] as a unit orthogonal vector, which is valid. So we can drop the requirement that \((x,y,z)\) must not be \((0,0,1)\). Moreover, we have \((s,t,u)=(0,0,1)\) only when \(z=0\), which is also a valid result. This leaves only the degeneracy from the denominator, i.e. \[ 1-x\cos\theta-y\sin\theta=0\,. \] (This is a removable discontinuity.) Since any rational formula for an orthogonal unit vector must agree with a formula presented above for some \(\theta\) almost everywhere, there is no formula that avoids at least one pesky edge case.Conclusion
We can see that Frisvad's formula is, up to renaming of variables, the expression we derived for \(\theta=0\) above. It is indeed the simplest, as the typical non-degenerate Pythagorean triples create expressions with significantly more terms. So unless chip manufacturers make some truly strange decisions about what they optimize on silicon, this formula is the way to go for orthogonal unit vectors.