This is a quick note about some methods for computing an orthonormal basis from a unit (bi)vector. The math directly follows that in “Implied normals (unit bivectors as rotations)”.


A re-cap of the needed bits from implied. Start with a standard orthonormal basis set:


Given a unit (bi)vector:


Then the quaternion $Q$ which rotates $\mathbf{z}$ to $\mathbf{v}$:


Toy C code can be found: HERE



Single reference direction sticking with up seems reasonable


Directly using $\eqref{fx}$ we can transform (rotate) the standard orthonormal set {$\mathbf{x,y,z}$} into a new set {$\mathbf{x’,y’,v}$}:


A simple translation of $\eqref{xp}$ and $\eqref{yp}$ into ‘C’ gives:


void ortho_basis_1(vec3_t* v, vec3_t* xp, vec3_t* yp)
{
  float x = -v->x;
  float y =  v->y;
  float z =  v->z; 
  float a = y/(z+1.f);  // y/(z+1)
  float b = y*a;        // y^2/(z+1)
  float c = x*a;        // -xy/(z+1)
  
  vec3_set(xp, z+b, c,      x);  // {z+y/(z+1),   -xy/(z+1), -x}
  vec3_set(yp, c,   1.f-b, -y);  // {-xy/(z+1), 1-y^2/(z+1), -y}
}


Simply inspecting the code or recalling the math shows a problem as $\mathbf{v}$ approaches $-\mathbf{z}$, we have an infinite number of solutions for $Q$ and numerically the math explodes. In cases were this is possible the above can be augmented by adding a branch and setting to xp and yp to two orthogonal vectors in the {$\mathbf{xy}$}-plane or the “correct” values in the limit. Choice of branch point and values depend on use requirements.


The math basis for this method when using positive $\mathbf{z}$ as a single reference direction is identically to that in Frisvad1. However the math is presented differently and the derivation is carried through differently as well:


which directly translation into code results at least four instead of two products.



Two opposite reference directions up and down, why not?


A second possible solution is to form the rotation with nearest direction of two choices. As an example we can form the rotation that maps $-\mathbf{z}$ to $\mathbf{v}$ and use that when $v_z < 0$:


Rotating $\mathbf{x}$ and $\mathbf{y}$ by $Q_n$:


The differences between $\eqref{xp} \eqref{yp}$ and $\eqref{xn} \eqref{yn}$ can be expressed by taking the absolute value of $v_z$ and products of $s_z= \text{sgn}\left(v_z\right)$. So we can combine the two expressions as follows:


static inline void ortho_basis_2(vec3_t* v, vec3_t* xp, vec3_t* yp)
{
  float x  = v->x;
  float y  = v->y;
  float z  = v->z;
  float sz = -sgn(z);
  float az = fabsf(z);
  float a  = y/(az+1.f); //   y/(|z|+1)
  float b  = y*a;        // y^2/(|z|+1)
  float c  = -x*a;       // -xy/(|z|+1)
  
  vec3_set(xp, az+b, c,     sz*x); // {|z|+y/(|z|+1),   -xy/(|z|+1), -x}
  vec3_set(yp, c,    1.f-b, sz*y); // {  -xy/(|z|+1), 1-y^2/(|z|+1), -y}
}


This method has a potential issue since $Q_n$ maps $-\mathbf{z}$ to $\mathbf{v}$. This means the source basis set is {$\mathbf{x,y,-z}$} when $v_z < 0$ which is left-handed2. This is not an issue if we simply need two mutually orthogonal normals in the complement plane, however for cases where it is an issue we can simply negate one of the computed normals in the $v_z < 0$ case which is equivalent to multiplying through by $s_z$. For both we then have:


Choosing to use $s_z \mathbf{x}’$ gives us:

void ortho_basis_2a(vec3_t* v, vec3_t* xp, vec3_t* yp)
{
  float x  = -v->x;
  float y  = v->y;
  float z  = v->z; 
  float sz = sgn(z);
  float a  = y/(fabsf(z)+1.f);
  float b  = y*a;
  float c  = x*a;
  
  vec3_set(xp, z+sz*b, sz*c, x); 
  vec3_set(yp, c,   1.f-b, -sz*y); 
}



Pixar’s method and inspiration thereof

addition added 20170401


A recently published paper from Pixar3 presents a new method directly derived from that of Frisvad $\eqref{frisvad_x} \eqref{frisvad_y}$ which is both efficient and high quality. Translating their branch-free version into equations we have:


Over reals these are equivalent to $\eqref{xpn}$ and $\eqref{ypnn}$ (so it returns the opposite results of the 2a method above). Their strategy for handling the negative half sphere is obviously superior to the one I used above. Where I was computing both abs and sign of z, they only use the latter. The transform is simple: multiply top and bottom through by $s_z$:


After the palmface I wanted to see what could be done with the new puzzle piece.

Starting with Pixar’s method we can negate both the returned vectors and on X64 hardware this lowers the number of issues from 36 to 324. Additionally both GCC and clang fail to remove the redundant sz*x term which bring us to 30:

void ortho_basis_pixar_r1(vec3_t* v, vec3_t* xp, vec3_t* yp)
{
  // -- snip start --
  float x  = v->x;
  float y  = v->y;
  float z  = v->z; 
  float sz = sgn(z);
  float a  = 1.0f/(sz+z);
  // -- snip end --
  float sx = sz*x;     // shouldn't be needed but is for gcc & clang, saves 2 issues
  float b  = x*y*a;
  vec3_set(xp, sx*x*a - 1.f, sz*b, sx);
  vec3_set(yp, b, y*y*a-sz, y);
}


For X86 we’ve already beaten method 2a in terms of issues, which requires 32.

If we can allow left-handed results on one of the half-spheres then we can multiply the $\mathbf{x}’$ result of the previous through by $s_z$ which lowers the number of X64 issues to 28 (same as lower quality method 2 above):

void ortho_basis_pixar_l1(vec3_t* v, vec3_t* xp, vec3_t* yp)
{
  // -- cut-n-paste snipped out part --
  float b  = x*y*a;
  vec3_set(xp, x*x*a - sz, b, x);
  vec3_set(yp, b, y*y*a - sz, y);
}


So far we haven’t changed the error. Since we are not allowing the compiler to treat FP as reals there are some redundant computations given the order of operations. Choosing to pull out y*a drops the number of issues to 29/26 respectively and my RMS measures5 show a minor negative impact6 of ~0.003%. (Both compilers figure out the sz*x terms here)

void ortho_basis_pixar_r2(vec3_t* v, vec3_t* xp, vec3_t* yp)
{
  // -- cut-n-paste snipped out part --
  float ya = y*a;
  float b  = x*ya;
  vec3_set(xp, sz*x*x*a - 1.f, sz*b, sz*x);
  vec3_set(yp, b, y*ya-sz, y);
}

void ortho_basis_pixar_l2(vec3_t* v, vec3_t* xp, vec3_t* yp)
{
  // -- cut-n-paste snipped out part --
  float ya = y*a;
  float b  = x*ya;
  vec3_set(xp, x*x*a - sz, b, x);
  vec3_set(yp, b, y*ya - sz, y);
}


That exhasts the things that seem profitable and jump out at me to try working with x component reduction used in $\eqref{frisvad_x}$.

Returning the my previous choice of x-component expansion: take $\eqref{xpnn} \eqref{ypnn}$, change to Pixar’s strategy and negating both give a left-handed system on one half sphere.

void ortho_basis_l1(vec3_t* v, vec3_t* xp, vec3_t* yp)
{
  float x  = v->x;
  float y  = v->y;
  float z  = v->z; 
  float sz = sgn(z);
  float a  = y/(z+sz);
  float b  = y*a;
  float c  = x*a;
  
  vec3_set(xp, -z-b, c, x);
  vec3_set(yp, c, b-sz, y);
}


and multiply $\mathbf{y}’$ by $s_z$ for right handed version:

void ortho_basis_r1(vec3_t* v, vec3_t* xp, vec3_t* yp)
{
  float x  = v->x;
  float y  = v->y;
  float z  = v->z; 
  float sz = sgn(z);
  float a  = y/(z+sz);
  float b  = y*a;
  float c  = x*a;
  vec3_set(xp, -z-b, c, x);
  vec3_set(yp, sz*c, sz*b-1, sz*y);
}


Note I’m only using X86 issues as a quick-n-dirty measure of goodness (easy to copy-n-past into godbolt).

A summary of results:

method x86-64 RMS rel-RMS notes
basis_2 28 2.5868×10-8 -7.028%  
basis_2a 32 2.5868×10-8 -7.028%  
pixar 36 2.4170×10-8 - as per paper
pixar_r1 30 2.4170×10-8 -  
pixar_r2 29 2.4215×10-8 -0.003%  
pixar_l1 28 2.4170×10-8 -  
pixar_l2 26 2.4215×10-8 -0.003%  
basis_r1 27 2.5868×10-8 -7.028%  
basis_l2 24 2.5868×10-8 -7.028%  



Other methods


Thanks to Stephen Hill for providing references to some existing methods and (better yet) for writting up a summary in a blog post7 and a link to a second pair of methods by Sam Hocevar8.

These methods are all roughly: Given $\mathbf{v}$ inspect the components and choose vector $\bar{\mathbf{w}}$ based on them. Cross product creates an orthgonal vector. This effectively partitions the sphere into “combing” groups.

The method above can be extended at increased computational complexity to shape how the normals are combed. As an example we could add a twist in {$\mathbf{x,y}$} like that found in a Wiechel projection or any of a number of craziness.



Viz

These figures roughly show the behavior of one output normal. I’m not finding a quick-and-dirty viz that I really like. Down positive z, positive x and negative z.

Method 1: Minimal angle orientation change on sphere. wtf

Method 2: Minimal angle orientation change on two half spheres. wtf

Hughes Moller: wtf


References and Footnotes

  1. “Building an orthonormal basis from a 3d unit vector without normalization”, Jeppe Frisvad, Journal of Graphics Tools 16:33, (2012) PDF/code page 

  2. If you don’t like “handedness” then $\mathbf{x}’ \times \mathbf{y}’ = -\mathbf{v}$ instead of $\mathbf{v}$. 

  3. “Building an Orthonormal Basis, Revisited”, Tom Duff, James Burgess, Per Christensen, Christophe Hery, Andrew Kensler, Max Liani, and Ryusuke Villemin, Journal of Computer Graphics Techniques Vol. 6, No. 1, 2017 JCGT 

  4. Instruction issue numbers are from x86-64 gcc 6.3 with -O3. x86-64 clang 4.0.0 -O3 numbers appear to be the same. 

  5. Take with a grain-of-salt. My RMS measures are not in agreement with those reported by Pixar. 

  6. Computed the standard way: $1-\frac{e}{e_0}$ 

  7. “Perpendicular Possibilities”, Stephen Hill, 2011 blog post 

  8. “On picking an orthogonal vector (and combing coconuts)”, Sam Hocevar, 2013 blog post