August 8th, 2017 Toy code can be found: HERE. All reported values are using clang without enabling floating-point contractions. Prelim the dead horse beating Skimming the math. Starting with a quaternion $Q=w+\left(x,y,z\right)$ then we can rotate $\mathbf{p}$ by: \[\begin{equation} \label{eq:sxform} \mathbf{p}' = Q\mathbf{p}Q^{-1} \end{equation}\] and if $Q$ is unit magnitude this reduces to: \[\begin{equation} \label{eq:rsxform} \mathbf{p}' = Q\mathbf{p}Q^* \end{equation}\] If $Q$ with magitude of $s$ is transformed by $\eqref{eq:rsxform}$ then the result is the composition of the rotation and uniform scaling of $s^2$. To create a matrix we need to apply the rotation to the basis set to form our three equations: \[\begin{eqnarray*} \mathbf{x} = \left(1,0,0\right) \\ \mathbf{y} = \left(0,1,0\right) \\ \mathbf{z} = \left(0,0,1\right) \\ \end{eqnarray*}\] which expanded and reduced gives: \[\begin{eqnarray} \mathbf{x}' & = & Q\mathbf{x}Q^{*} \nonumber \\ & = & \left(w^2 + x^2 - y^2 - z^2, ~2\left(xy+wz\right), ~2\left(xz-wy\right) \right) \label{xp} \\ & = & \left(1 - 2 \left(y^2+z^2\right), ~2\left(xy+wz\right), ~2\left(xz-wy\right) \right) \label{xpr} \\ \nonumber \\ \mathbf{y}' & = & Q\mathbf{y}Q^{*} \nonumber \\ & = & \left(2 \left(xy-wz\right), ~w^2 - x^2 + y^2 - z^2, ~2\left(wx+yz\right) \right) \label{yp} \\ & = & \left(2 \left(xy-wz\right), ~1 - 2\left(x^2+z^2\right), ~2\left(wx+yz\right) \right) \label{ypr} \\ \nonumber \\ \mathbf{z}' & = & Q\mathbf{z}Q^{*} \nonumber \\ & = & \left(2 \left(wy+xz\right), ~2\left(yz-wx\right), ~w^2 - x^2 - y^2 + z^2\right) \label{zp} \\ & = & \left(2 \left(wy+xz\right), ~2\left(yz-wx\right), ~1 - 2\left(x^2+y^2\right) \right) \label{zpr} \end{eqnarray}\] Quaternion to rotation matrix Sticking to the math convention of column vectors, then we can shove the (not reduced) equations $\eqref{xp} \eqref{yp} \eqref{zp}$ into the columns giving: \[\begin{equation} \label{2matndr} { \left( \begin{array}{ccc} w^2 + x^2 - y^2 - z^2 & 2\left(xy-wz\right) & 2\left(xz+wy\right) \\ 2 \left(xy+wz\right) & w^2 - x^2 + y^2 - z^2 & 2\left(yz-wx\right) \\ 2 \left(xz-wy\right) & 2\left(yz+wx\right) & w^2 - x^2 - y^2 + z^2 \end{array} \right) } \end{equation}\] and we can reduce the diagonal elements since $Q\cdot Q = 1$ as per $\eqref{xpr} \eqref{ypr} \eqref{zpr}$: \[\begin{equation} \label{2matstd} { \left( \begin{array}{ccc} 1 - 2 \left(y^2+z^2\right) & 2\left(xy-wz\right) & 2\left(xz+wy\right) \\ 2 \left(xy+wz\right) & 1 - 2\left(x^2+z^2\right) & 2\left(yz-wx\right) \\ 2 \left(xz-wy\right) & 2\left(yz+wx\right) & 1 - 2\left(x^2+y^2\right) \end{array} \right) } \end{equation}\] Translating both into C gives: // not reducing the diagonal. adds ~4 (basic simple scalar) ops // vs. the standard version. void quat_to_mat33_ndr(mat33_t* m, quat_t* q) { float x = q->x, y = q->y, z = q->z, w = q->w; float tx = 2*x, ty = 2*y, tz = 2*z; float xy = ty*x, xz = tz*x, yz = ty*z; float wx = tx*w, wy = ty*w, wz = tz*w; // diagonal terms float t0 = (w+y)*(w-y), t1 = (x+z)*(x-z); float t2 = (w+x)*(w-x), t3 = (y+z)*(y-z); m->m00 = t0+t1; m->m11 = t2+t3; m->m22 = t2-t3; m->m10 = xy+wz; m->m01 = xy-wz; m->m20 = xz-wy; m->m02 = xz+wy; m->m21 = yz+wx; m->m12 = yz-wx; } // commonly seen variant void quat_to_mat33_std(mat33_t* m, quat_t* q) { float x = q->x, y = q->y, z = q->z, w = q->w; float tx = 2 *x, ty = 2 *y, tz = 2 *z; float xx = tx*x, yy = ty*y, zz = tz*z; float xy = ty*x, xz = tz*x, yz = ty*z; float wx = tx*w, wy = ty*w, wz = tz*w; m->m00 = 1.f-(yy+zz); m->m11 = 1.f-(xx+zz); m->m22 = 1.f-(xx+yy); m->m10 = xy+wz; m->m01 = xy-wz; m->m20 = xz-wy; m->m02 = xz+wy; m->m21 = yz+wx; m->m12 = yz-wx; } NOTE: An identical problem is extracting one or more axes from a quaterion. Say $\mathbf{x}, \mathbf{y}, \mathbf{z}$ respectively represent: right, forward, up. We can find the transformed direction(s) using the corresponding primed equation(s). The reason for this aside that because I've sometime seen code that extracts two of the directions as above and builds the third using a cross product. For the cross product method to be a performance win requires that it ends up being faster than three additions. Another aside: In cases where input quaternions are moving away from unity due to compounding errors it is possible to fold the corrective step into the transform instead of renormalizing as a pre-step. Rotation matrix to quaternion Let’s rename the elements of $\eqref{2matndr}$ as: \[{ \left( \begin{array}{ccc} m_{00} & m_{01} & m_{02} \\ m_{10} & m_{11} & m_{12} \\ m_{20} & m_{21} & m_{22} \end{array} \right) }\] then the sum and difference of off diagonal pairs give: \[\begin{eqnarray} m_{10} + m_{01} & = & 4xy \label{4xy} \\ m_{10} - m_{01} & = & 4wz \label{4wz} \\ m_{02} + m_{20} & = & 4xz \label{4xz} \\ m_{02} - m_{20} & = & 4wy \label{4wy} \\ m_{21} + m_{12} & = & 4yz \label{4yz} \\ m_{21} - m_{12} & = & 4wx \label{4wx} \end{eqnarray}\] and the permutations of trace like functions: \[\begin{eqnarray*} m_{00} + m_{11} + m_{22} & = & 3 w^2 - \left(x^2 + y^2 + z^2 \right) & = & -1+4w^2 \\ m_{00} + m_{11} - m_{22} & = & \left(w^2 + x^2 + z^2 \right)- 3 y^2 & = & 1-4y^2 \\ m_{00} - m_{11} + m_{22} & = & \left(w^2 + x^2 + y^2 \right) - 3 z^2 & = & 1-4z^2 \\ m_{00} - m_{11} - m_{22} & = & 3 x^2 - \left(w^2 + y^2 + z^2\right) & = & -1+4x^2 \\ \end{eqnarray*}\] which can be reworked to: \[\begin{eqnarray} 4w^2 & = & 1 + m_{00} + m_{11} + m_{22} \label{4w2} \\ 4x^2 & = & 1 + m_{00} - m_{11} - m_{22} \label{4x2}\\ 4y^2 & = & 1 - m_{00} + m_{11} - m_{22} \label{4y2}\\ 4z^2 & = & 1 - m_{00} - m_{11} + m_{22} \label{4z2} \end{eqnarray}\] If we know or determine an element has a sufficiently large magnitude we can choose the equation for that element from the second set of equations and three equations from the first set which are products of the that and the remaining elements. As an example take a special case version where we know that the magnitude of the angle is not greater than $\frac{\pi}{2}$ which in turn ensures that the scalar $(w)$ is the largest element. Then we need $\eqref{4w2}$ from the second set and $\eqref{4wx} \eqref{4wy} \eqref{4wz}$ from the first. void mat33_2_quat_small(quat_t* q, mat33_t* m) { float m00=m->m00, m01=m->m01, m02=m->m02; float m10=m->m10, m11=m->m11, m12=m->m12; float m20=m->m20, m21=m->m21, m22=m->m22; float t = 1 + m00 + m11 + m22; // 4w^2 float s = 0.5f*rsqrtf(t); // 1/(4w) quat_set(q, s*(m21-m12), s*(m02-m20), s*(m10-m01), s*t); } Standard construction The standard construction for the general case is based on the observation that since we’re producing a unit quaternion, then the smallest that the largest magnitude component can be is $\frac{1}{2}$. Given that we can simply walk through the components and choose the first that’s at least that large. void mat33_to_quat_std(quat_t* q, mat33_t* m) { float m00=m->m00, m01=m->m01, m02=m->m02; float m10=m->m10, m11=m->m11, m12=m->m12; float m20=m->m20, m21=m->m21, m22=m->m22; float t1 = m22+m11; float t0 = m00+t1; if (t0 > 0.f) { float w = 1.f+t0, s = 0.5f*rsqrtf(w); quat_set(q, s*(m21-m12), s*(m02-m20), s*(m10-m01), s*w); // w return; } t0 = m00-t1; if (t0 > 0.f) { float x = 1.f+t0, s = 0.5f*rsqrtf(x); quat_set(q, s*x, s*(m01+m10), s*(m02+m20), s*(m21-m12)); // x return; } t0 = m11-m22; t1 = 1.f-m00; if (t0 > 0.f) { float y = t1+t0, s = 0.5f*rsqrtf(y); quat_set(q, s*(m01+m10), s*y, s*(m12+m21), s*(m02-m20)); // y } else { float z = t1-t0, s = 0.5f*rsqrtf(z); quat_set(q, s*(m02+m20), s*(m12+m21), s*z, s*(m10-m01)); // z } } Given uniform random input then the (in code order) probabilities are: $\frac{2}{3},\frac{1}{9},\frac{1}{9},\frac{1}{9} $ The first from $\frac{2}{\pi}\text{acos}\left(\frac{1}{2}\right) = \frac{2}{3}$. Mike Day’s construction Insomniac Games (Mike Day) presented a paper1 in 2015 which makes the following observations. The sign of diagonal components imply $\eqref{2matstd}$: \[\begin{eqnarray*} m_{00} < 0 \implies y^2+z^2 > \frac{1}{2} \\ m_{11} < 0 \implies x^2+z^2 > \frac{1}{2} \\ m_{22} < 0 \implies x^2+y^2 > \frac{1}{2} \end{eqnarray*}\] and the permutations of sum/difference of pairs of diagonal elements imply: \[\begin{eqnarray*} m_{00}+m_{11} < 0 \implies w^2 < z^2 \\ m_{00}-m_{11} < 0 \implies x^2 < y^2 \\ m_{00}+m_{22} < 0 \implies w^2 < y^2 \\ m_{00}-m_{22} < 0 \implies x^2 < z^2 \\ m_{11}+m_{22} < 0 \implies w^2 < x^2 \\ m_{11}-m_{22} < 0 \implies y^2 < z^2 \end{eqnarray*}\] They use the above to form a method that insures to choose a sufficiently large element. The following is a slightly massaged version of Day: void mat33_2_quat_day(quat_t* q, mat33_t* m) { float m00=m->m00, m01=m->m01, m02=m->m02; float m10=m->m10, m11=m->m11, m12=m->m12; float m20=m->m20, m21=m->m21, m22=m->m22; float e0; if (m22 >= 0) { float a = m00+m11; float b = m10-m01; float c = 1.f+m22; if (a >= 0) { e0 = c + a; quat_set(q, m21-m12, m02-m20, b, e0); // w } else { e0 = c - a; quat_set(q, m02+m20, m21+m12, e0, b); // z } } else { float a = m00-m11; float b = m10+m01; float c = 1.f-m22; if (a >= 0) { e0 = c + a; quat_set(q, e0, b, m02+m20, m21-m12); // x } else { e0 = c - a; quat_set(q, b, e0, m21+m12, m02-m20); // y } } quat_scale(q, 0.5f*rsqrtf(e0)); } Given uniform random input then the outer and both inner probabilities are: $\frac{1}{2}$. Angle errors Let’s examine round-trip errors. We generate a uniform quaternion $A$, convert to a matrix, back to a quaternion $B$ and then find the relative quaternion $X$: \[X = BA^* = \mathbf{v} + w\] then measure the represented angle, multiple by two (to get 3D rotation measure instead of in $\mathbb{H}$) and convert to degrees: \[\frac{360}{\pi}\text{atan}\left(\frac{\left\Vert \mathbf{v} \right\Vert }{\abs{w}}\right)\] This choice of measure ignores the magnitudes of $A$ and $B$. The trace names are to-matrix-method/to-quat-method. The x-axis is normalized. There’s a jump in trace behavior at $\frac{2}{3}$ because our x-axis is a function of $w$ and this the point where $w=\frac{1}{2}$. The forward transform without diagonal reduction has lower peak error, but potentially more interesting is the fact that this method is better at angle-preserving (conformal) than the standard. Variants you probably don’t care Branch free Day The reworking of Day above is to point out the we can formulate a branch-free version. Let’s call a function to isolate the sign bit of a float $sb$ and denote an integer XOR as $\oplus$. The initial branch is selected based on the sign of $m_{22}$. In both cases we’re computing three temp variables $a,b,c$ which can be expressed in terms of that sign: $s_0 = \text{sb}\left(m_{22}\right)$. Then inside the inner four cases we’re computing $e_0$ (the chosen sufficiently large element) which can be expressed in terms of the sign of $a$: $s_1 = \text{sb}\left(a\right)$. ($c$ and $e_0$ can use abs instead). The remaining two elements (unnamed temp expressions) can be formed using $s_1$ and the parity of two sign values: $s_p = s_0 \oplus s_1$. If we rename $b$ to $e_1$, then so far we have: \[\begin{eqnarray*} s_0 & = & \text{sb}\left(m_{22}\right) \\ a & = & m_{00} + \left(m_{11} \oplus s_0\right) \\ c & = & 1 + \left(m_{22} \oplus s_0\right) = 1 + \abs{m_{22}} \\ s_1 & = & \text{sb}\left(a\right) \\ s_p & = & s_0 \oplus s_1 \\ e_0 & = & c + \left(a \oplus s_1\right) = c + \abs{a} \\ e_1 & = & m_{10} - \left(m_{01} \oplus s_0\right) \\ e_2 & = & m_{02} - \left(m_{20} \oplus s_p\right) \\ e_3 & = & m_{21} - \left(m_{12} \oplus s_1\right) \end{eqnarray*}\] To complete the component computation we simply need to scale all four by $\frac{1}{2 \sqrt{e_0}}$. This leave placing them in memory, if we has storage order of $\left(x,~y,~z,~w\right)$ and think of our elements as the ordered list: $\left(e_3,~e_2,~e_1,~e_0 \right) $, then we apply the follow two rules (input to the second is renaming the output of the first) if $s_0$ is set $\left(e_3,~e_2,~e_1,~e_0 \right) \rightarrow \left( e_1,~e_0,~e_3,~e_2 \right) $, so flip the first and last two. if $s_p$ is set $\left(e_3,~e_2,~e_1,~e_0 \right) \rightarrow \left( e_2,~e_3,~e_0,~e_1 \right) $, so flip the first and second pairs. Alternately we can simply compute a value to XOR as below: void mat33_to_quat_day_bf(quat_t* q, mat33_t* m) { float m00=m->m00, m01=m->m01, m02=m->m02; float m10=m->m10, m11=m->m11, m12=m->m12; float m20=m->m20, m21=m->m21, m22=m->m22; // isolate a sufficently large component uint32_t s0 = sb(m22); float a = m00 + fxor(m11,s0); // m00 + (m11^s0) uint32_t s1 = sb(a); uint32_t sp = s0^s1; float e0 = 1.f + fxor(m22,s0) + fxor(a,s1); // 1 + |m22| + |a| float s = 0.5f*rsqrtf(e0); // the other three components float e1 = m10-fxor(m01,s0); // m10 - (m01^s0) float e2 = m02-fxor(m20,sp); // m02 - (m02^sp) float e3 = m21-fxor(m12,s1); // m21 - (m12^s1) e0 *= s; e1 *= s; e2 *= s; e3 *= s; // scaling uint32_t id = (s0 >> 30)^(sp >> 31); // write ordering q->f[3^id] = e0; q->f[2^id] = e1; q->f[1^id] = e2; q->f[0^id] = e3; } Branch free higher error Another branch-free possibility is to use the second set of equations to compute the magnitudes of all the elements, choose the scalar to be positive and compute the sign of the bivector components by the first set of equations. void mat33_2_quat_bf0(quat_t* q, mat33_t* m) { float m00=m->m00, m01=m->m01, m02=m->m02; float m10=m->m10, m11=m->m11, m12=m->m12; float m20=m->m20, m21=m->m21, m22=m->m22; float a = 1.f+m00; float b = 1.f-m00; float c = m11+m22; float d = m11-m22; // abs to prevent errors, could use max against zero // or a bias amount as per below. float w = fabsf(a+c); // 4w^2 float x = fabsf(a-c); // 4x^2 float y = fabsf(b+d); // 4y^2 float z = fabsf(b-d); // 4z^2 // assumes rsqrtf(0) = inf, otherwise needs to add a bias w *= 0.5f*rsqrtf(w); x *= 0.5f*rsqrtf(x); y *= 0.5f*rsqrtf(y); z *= 0.5f*rsqrtf(z); // use first set for the signs x = copysignf(x, m21-m12); y = copysignf(y, m02-m20); z = copysignf(z, m10-m01); quat_set(q, x,y,z,w); } which as promised results in much higher errors: Although there are cases where these errors are low enough to be usable the performace is worse that some better methods… Almost branch free Yet another option is to notice that special case small angle version performs rather well. These plots again are: So we can formulate a version that has a highly predictable branch that handles the common case and otherwise fall-back. void mat33_to_quat_bf1(quat_t* q, mat33_t* m) { float m00=m->m00, m01=m->m01, m02=m->m02; float m10=m->m10, m11=m->m11, m12=m->m12; float m20=m->m20, m21=m->m21, m22=m->m22; float t = 1.f + m00 + m11 + m22; // 4w^2 if (t > BF1_CUT) { float s = 0.5f*rsqrtf(t); quat_set(q, s*(m21-m12), s*(m02-m20), s*(m10-m01), s*t); return; } // add handling other cases here } A lazy table of empirical numbers: some ad-hoc cut points, probability and round-trip errors using the standard method of quat-to-matrix: BF1_CUT probability rt-error 0.00001 0.998002 0.042664 0.0001 0.993634 0.014877 0.0002 0.990994 0.010233 0.0004 0.987265 0.007769 0.0008 0.981993 0.004991 0.0016 0.974531 0.003587 0.0032 0.963969 0.002935 0.0064 0.949045 0.001946 0.0128 0.927966 0.001367 0.0256 0.898184 0.000975 0.0512 0.856216 0.000753 0.1024 0.797129 0.000528 Overkill but lower error Let’s run through a method which includes every component of the rotation matrix in each component of the resulting quaternion. The point of that exercise is a light filtering of noise present in the input. First we’ll rewrite the rotation matrix as the result of converting an axis-angle conversion with unit vector $(x,~y,~z)$ and angle $\theta \in [0,\pi]$: \[\left( \begin{array}{ccc} \cos \theta +x^2 (1-\cos \theta ) & x y (1-\cos \theta )-z \sin \theta & x z (1-\cos \theta)+y \sin \theta \\ x y (1-\cos \theta )+z \sin \theta & \cos \theta +y^2 (1-\cos \theta ) & y z (1-\cos \theta )-x \sin \theta \\ x z (1-\cos \theta )-y \sin \theta & y z (1-\cos \theta) + x \sin \theta & \cos \theta +z^2 (1-\cos \theta ) \\ \end{array} \right)\] NOTE: everywhere else I’m using $x,y,z$ for quaternion components and here for the unit vector of the axis angle. Expanding and naming off diagonal differences: \[\begin{eqnarray*} s_1 = m_{21} - m_{12} = 2x \sin(\theta) \\ s_2 = m_{02} - m_{20} = 2y \sin(\theta) \\ s_3 = m_{01} - m_{10} = 2z \sin(\theta) \end{eqnarray*}\] and sums: \[\begin{eqnarray*} a_1 = m_{21} + m_{12} = 2yz \left(1-\cos \theta \right) \\ a_2 = m_{02} + m_{20} = 2xz \left(1-\cos \theta \right) \\ a_3 = m_{01} + m_{10} = 2xy \left(1-\cos \theta \right) \end{eqnarray*}\] and (what I’m calling) trace like: \[\begin{array}{lll} t_0 = & m_{00} + m_{11} + m_{22} = & 1 + 2~\cos \theta \\ t_1 = & m_{00} - m_{11} - m_{22} = & 2~x^2 \left(1-\cos \theta \right)-1 \\ t_2 = & m_{11} - m_{00} - m_{22} = & 2~y^2 \left(1-\cos \theta \right)-1 \\ t_2 = & m_{22} - m_{00} - m_{11} = & 2~z^2 \left(1-\cos \theta \right)-1 \end{array}\] We can combine these as follows: \[\begin{array}{lll} r_0 = (t_0+1)^2 + s_1^2 + s_2^2 + s_3^2 = & 8\left(1+\cos \theta \right) \\ r_1 = (t_1+1)^2 + s_1^2 + a_2^2 + a_3^2 = & 8\left(x^2 \left(1-\cos \theta \right)\right) \\ r_2 = (t_2+1)^2 + s_2^2 + a_1^2 + a_3^2 = & 8\left(y^2 \left(1-\cos \theta \right)\right) \\ r_3 = (t_3+1)^2 + s_3^2 + a_1^2 + a_2^2 = & 8\left(z^2 \left(1-\cos \theta \right)\right) \end{array}\] a simple rewrite to make the half-angle identities jump out: \[\begin{array}{lll} r_0 = 16 \left(\frac{1+\cos \theta}{2} \right) \\ r_1 = 16~x^2 \left(\frac{1-\cos \theta}{2} \right) \\ r_2 = 16~y^2 \left(\frac{1-\cos \theta}{2} \right) \\ r_3 = 16~z^2 \left(\frac{1-\cos \theta}{2} \right) \end{array}\] So we can compute the magnitudes of the quaternion components as follows: \[\begin{array}{lll} |q_w| = \frac{\sqrt{r_0}}{4} \\ |q_x| = \frac{\sqrt{r_1}}{4} \\ |q_y| = \frac{\sqrt{r_2}}{4} \\ |q_z| = \frac{\sqrt{r_3}}{4} \end{array}\] The remaining issue is then determining the signs of the components and the $s_i$ terms have that information. Plot the error of that implementation (vs. Day) would give us: So it all falls apart when approaching the maximal rotation angle. The problem is $\sin \theta$ is approaching zero. We can correct this by noting the most important information is in the largest magnitude component $q_i$ and we could choose that to be the one that’s positive (or equivalently known sign). Specifically the $a_i$ term which includes the known and desired sign components. With this correction we have: Let’s do some points: dp is promoting input to doubles, computing and demoting to single. For single precision this is less work than being fancy. If working precision is doubles you’d have to go with extended precision methods (TwoProduct, compensated summations, etc.) to squeeze out the additonal error reduction. Although it’s structure helps with noise filtering, performing the standard quaternion to matrix method (reducing the diagaonals) is still a major poke in the eye. Just hide all the std methods now..thanks! All of this added computational complexity does buy us anything for small rotation angles (vs. Day & “Almost branchfree”..provided NDR quat-to-mat) Since NDR is the most effective error reduction, let’s do a spitball of it promoted to double computation. (reminder: the toy code will let you play with whatever combo you want) Now we’re not really attempting to filter a noisy input we can relax to ensuring $w$ is sufficiently large to not increase the error in the case we’re interested it to not have a branch misprediction nightmare. (Going branchfree is possible and would be the “thing” todo in a SIMD implementation. For scalar it doesn’t seem worthwhile). Making the simplifing assumption that rotations are uniformly random distributed then the normalized cumulative volume (cumulative/total) of a rotation is: \[\frac{\theta - \sin \theta}{\pi}\] for a 3D rotation angle of $\theta$. Plugging in how the quaternion scalar $w$ is related to $\theta = 2 ~\text{acos} (w)$ yields: \[\frac{2 ~\text{acos} (w) - 2w\sqrt{1-w^2}}{\pi}\] which is pretty close to a linear decreasing function: So if the first test is w > CUT then the probablity of that block is $w$ fed into the above expression (assuming uniform input and better if angles are skewed toward small). An implementation: // example cut points: // @ 0.1015625 reaches 'w' case 87.1% <- this is playin' it safe // @ 0.0678711 reaches 'w' case 91.3% <- this is close to min // value w/o error increase // example choice #define OKC_CUT 0.0678711f void mat33_to_quat_okc(quat_t* q, mat33_t* m) { float m00=m->m00, m01=m->m01, m02=m->m02; float m10=m->m10, m11=m->m11, m12=m->m12; float m20=m->m20, m21=m->m21, m22=m->m22; float t0 = m00+m11+m22, d0 = t0+1.f; float t1 = m00-m11-m22, d1 = t1+1.f; float t2 = m11-m00-m22, d2 = t2+1.f; float t3 = m22-m00-m11, d3 = t3+1.f; float s1 = m21-m12, s2 = m02-m20, s3 = m10-m01; float a1 = m21+m12, a2 = m02+m20, a3 = m10+m01; float w = 0.25f*sqrtf(d0*d0 + s1*s1 + s2*s2 + s3*s3); float x = 0.25f*sqrtf(d1*d1 + s1*s1 + a2*a2 + a3*a3); float y = 0.25f*sqrtf(d2*d2 + s2*s2 + a1*a1 + a3*a3); float z = 0.25f*sqrtf(d3*d3 + s3*s3 + a1*a1 + a2*a2); // could be make branchfree (thinkin' SIMD here) if (w > OKC_CUT) { x = copysignf(x, s1); y = copysignf(y, s2); z = copysignf(z, s3); } else { if (x > 0.25f) { y = copysignf(y,a3); z = copysignf(z,a2); w = copysignf(w,s1); } else if (y > 0.25f) { x = copysignf(x, a3); z = copysignf(z, a1); w = copysignf(w, s2); } else { x = copysignf(x, a2); y = copysignf(y, a1); w = copysignf(w, s3); } } quat_set(q, (float)x,(float)y,(float)z,(float)w); } References and Footnotes “Converting a Rotation Matrix to a Quaternion”, Mike Day, 2015 (download page) ↩ Comments quaternions (15) rotations (9) , conversion (1)