How on Earth one converts 3x3 matrix to quaternion?

Unity Engine
,
1

Edit: This was solved. Quaternion.LookRotation C# (re-)implementation (+ standalone 3x3 orthogonal rotation matrix → quaternion conversion) can be found in post #9 . + bonus content: Quaternion.ToAngleAxis in post #12 .

Heya,

First of all, I’m doing this for learning purposes. Please don’t instruct me to not reinvent the wheel.

Next, before jumping to conclusions, let me say this: I’ve tried everything one can find on Google. Or ChatGPT. Or Copilot. Or Wikipedia. Or Siggraph papers. Coding gems. Decompiled Unity source. You name it. It is all wrong. Or I’m not getting something big time.

Of course, it’s extremely likely that I’m missing a step here and tbh I’m not that versed in eigenvalues and matrix decomposition, SO(3), hyperspheres, Lie algebras etc. Not to mention the handedness jungle, randomly transposed matrices, various quaternion standards (ijkw, wijk, wxyz, xyzw …), as well as Unity itself being non-transparent about it…

I’m honestly a dum-dum when it comes to “ah you just do #&^# @&^@# and then because w is octohugging the perplexdome of the verticultural jynxopsis you can #$&# @&# and look how easy is that, you arrive at ∫Z(x^2) e^ipi and that’s your solution.”

I’m sorry but if I wanted to feel depressed I’d play Starfield instead.

Anyway I also checked out the source for Unity mathematics, but this particular thing (in float3x3) looks like pizza from hell. It’s impossible to tell what’s going on with all the conversions from uint to float, extracting signs then XORing them back…

These two are correct however, so I have something to lean onto

var r = Quaternion.LookRotation(fwd, up);

and

Vector3 rgt = cross(up, fwd);
up = cross(fwd, rgt);
var r = new Matrix4x4(rgt, up, fwd, new Vector4(0f, 0f, 0f, 1f)).rotation;

(Both of these will jump to the C++ side for the actual conversion and I believe it’s the exact same code.)

And so obviously I used the latter as a basis for my analysis, because the actual matrix is clear as day.

| Rx Ux Fx 0 |
| Ry Uy Fy 0 |
| Rz Uz Fz 0 |
|  0  0  0 1 |

Here’s what I tried:
If I learned anything it’s that the quaternion is obtained by extracting the ‘trace of the matrix’ T

T = Rx + Uy + Fz;

And then w = sqrt(1 + T) / 2

Most of the implementations I’ve seen have branching where they take some piece of the matrix with the “biggest meat on it” and reshuffle all components in order to avoid divisions by zero. I’ve also seen many odd takes that claim to be numerically stable etc, on StackOverflow and in Unity Answers, but when you take it all in, some are just naive, or somehow messed up, and many are literal rewrites of the same thing over and over.

After some time, I’ve grown impatient with trying to find the whole thing in a pristine condition, so I’ve decided to probe the values themselves in a contrived scenario.

First, let’s say we don’t care about the signs. Then we note that every single implementation uses the above formula for exactly one of the final components (of course T is true only for w and the expression will change for x, y, z, so let’s call this N instead).

Because this evaluation only needs a diagonal, it doesn’t matter whether I’ve transposed the matrix (so at least that’s out). And it should be guaranteed that the value of sqrt(1 + N) / 2 must appear somewhere in the quaternion. However there are 4 possible expressions for N

N =  Rx + Uy + Fz; // w (T)
N =  Rx - Uy - Fz; // x
N = -Rx + Uy - Fz; // y
N = -Rx - Uy + Fz; // z

This is the scenario I came up with:
With some Target at (2, 8, 2), we get the following “LookRotation” basis vectors

R ( 0.23570, 0.00000, -0.23570)
U (-0.22222, 0.11111, -0.22222)
F ( 0.23570, 0.94281,  0.23570)

How to arrive here

var target = new Vector3(2f, 8f, 2f);
var fwd = target.normalized; // unit
var up = Vector3.up; // unit
var rgt = cross(up, fwd); // Edit: I should've normalized this; it's so funny I didn't write 'unit' here -- my brain apparently knew it all along
up = cross(fwd, rgt); // re-orthogonalized
Debug.Log($"R {rgt:F5}");
Debug.Log($"U {up:F5}");
Debug.Log($"F {fwd:F5}");

The diagonal of this matrix (Rx, Uy, Fz) (or in another notation m00, m11, m22) is thus
(0.23570, 0.11111, 0.23570)

LookRotation (and the other Matrix4x4.rotation example) will return this quaternion
(-0.53340, 0.31246, 0.22094, 0.75434)
How to arrive here

var q = Quaternion.LookRotation(new Vector3(2f, 8f, 2f).normalized, Vector3.up);
Debug.Log(q);

But when I ran sqrt(1 + N) / 2 (for all N as shown above) none of the values I got from it matched with the bolded values. NONE! To be safe, I included 4 more sign permutations (for a complete set of 2^3=8) and nothing. I also tried normalizing the quaternion, which shouldn’t be required, to no avail. I tried as well permuting the signs on output or changing sqrt(1 + N) / 2 to absurdity. I don’t even know why I tried that, I guess I’m at the point of randomizing all of the things I tried so far in hopes of stumbling into some miraculous insight.

So the question is, how on Earth one obtains the bolded values above?

Apart from the obvious results in Google, here’s a useful PDF that mentions two of the most prominent solutions I’ve seen before (and any AI or Quora bots will readily recite these like a psalm). I’ve tried numerous variations of the both examples shown here and nothing made any sense. Edit: there is also this on Wikipedia, that states the exact same thing.

Where is my error?

Edit: Fixed trace of matrix from 1+Rx+Uy+Fz to just Rx+Uy+Fz as it should be (1 is pulled out to sqrt(1 + N) / 2). Nothing changes because of it, I just thought it would look more correct.

Edit2: After I found the bug, I added some remarks above. The issue was that I should’ve normalized the right vector because up and forward aren’t guaranteed to be orthogonal at that moment.

2

I don’t know how you get the bolded values above, sorry. But to get a quaternion from a 4x4 matrix containing a valid rotation matrix I think you can just use:

You can form a 4x4 from a 3x3 as you have shown above. 3x3 rotation matrices have special conditions on them, not all 3x3 matrices are valid rotations. Good luck!

1 Like
3

I don’t think orionsyndrome is asking this in terms of a function in the API.

But rather is trying to learn the actual logic behind the conversion.

…

Sorry orion, I’d love to hang around and assist right now. But with the holiday my sched can’t afford the time. If I find a break I may come back… or just to read what other’s have to say.

2 Likes
4

Yes, eigen decomposition is something quite complex. It’s sad that 3b1b’s Linear algebra series only mentions what eigenvalues / eigenvectors are, but not how to calculate them :slight_smile: The Wikipedia article mentioned that a 3x3 rotation matrix has at least one eigenvalue with the value 1, the others are actually complex numbers (because the trace / square root is negative). The article also used the example of using the usual trace of the matrix and just states that if it’s zero it doesn’t work -.-

The actual vector is essentially similar to the crossproduct. So the x component is always the result of y and z, y is the result of x and z and finally z is the result of x and y. Though choosing which components seems to be a bit tricky. I found this post on Unity Answers / Discussions which seems to have a working version of LookRotation (which starts by constructing the 3x3 matrix).

I took your example direction and made this test case with the method I found on Discussions:

Test Example

using UnityEngine;

public class QuaternionTest : MonoBehaviour
{
    void Start()
    {
        var target = new Vector3(2f, 8f, 2f);
        var fwd = target.normalized;
        var up = Vector3.up;
        var rgt = Vector3.Cross(up, fwd);
        up = Vector3.Cross(fwd, rgt);

        var q1 = Quaternion.LookRotation(fwd, up);
        Debug.Log(q1.ToString("f8"));

        var q2 = QuaternionLookRotation(fwd, up);
        Debug.Log(q2.ToString("f8"));

    }

    private static Quaternion QuaternionLookRotation(Vector3 forward, Vector3 up)
    {
        forward.Normalize();

        Vector3 vector = Vector3.Normalize(forward);
        Vector3 vector2 = Vector3.Normalize(Vector3.Cross(up, vector));
        Vector3 vector3 = Vector3.Cross(vector, vector2);
        var m00 = vector2.x;
        var m01 = vector2.y;
        var m02 = vector2.z;
        var m10 = vector3.x;
        var m11 = vector3.y;
        var m12 = vector3.z;
        var m20 = vector.x;
        var m21 = vector.y;
        var m22 = vector.z;


        float num8 = (m00 + m11) + m22;
        var quaternion = new Quaternion();
        if (num8 > 0f)
        {
            var num = (float)System.Math.Sqrt(num8 + 1f);
            quaternion.w = num * 0.5f;
            num = 0.5f / num;
            quaternion.x = (m12 - m21) * num;
            quaternion.y = (m20 - m02) * num;
            quaternion.z = (m01 - m10) * num;
            return quaternion;
        }
        if ((m00 >= m11) && (m00 >= m22))
        {
            var num7 = (float)System.Math.Sqrt(((1f + m00) - m11) - m22);
            var num4 = 0.5f / num7;
            quaternion.x = 0.5f * num7;
            quaternion.y = (m01 + m10) * num4;
            quaternion.z = (m02 + m20) * num4;
            quaternion.w = (m12 - m21) * num4;
            return quaternion;
        }
        if (m11 > m22)
        {
            var num6 = (float)System.Math.Sqrt(((1f + m11) - m00) - m22);
            var num3 = 0.5f / num6;
            quaternion.x = (m10 + m01) * num3;
            quaternion.y = 0.5f * num6;
            quaternion.z = (m21 + m12) * num3;
            quaternion.w = (m20 - m02) * num3;
            return quaternion;
        }
        var num5 = (float)System.Math.Sqrt(((1f + m22) - m00) - m11);
        var num2 = 0.5f / num5;
        quaternion.x = (m20 + m02) * num2;
        quaternion.y = (m21 + m12) * num2;
        quaternion.z = 0.5f * num5;
        quaternion.w = (m01 - m10) * num2;
        return quaternion;
    }
}

This is the result:
9922137--1435146--upload_2024-7-3_22-33-53.png

So it seems to generate the same result. Now we should test this with random directions to see if it actually works out in all (most) cases.

1 Like
5

I just changed the Start method of my test code to this:

    void Start()
    {
        int success = 0;
        int max = 1000;
        for (int i = 0; i < max; i++)
        {
            var target = Random.onUnitSphere;
            var fwd = target.normalized;
            var up = Random.onUnitSphere;
            var rgt = Vector3.Cross(up, fwd);
            up = Vector3.Cross(fwd, rgt);

            var q1 = Quaternion.LookRotation(fwd, up);
            var q2 = QuaternionLookRotation(fwd, up);
            if (q1 != q2)
                Debug.Log($"Not equal {q1.ToString("f8")} - {q2.ToString("f8")}");
            else
                success++;
        }
        Debug.Log($"{success} / {max}");
    }

Here I check 1000 completely random rotation matrices and I get a success rate of 100%. So not one case where q1 and q2 do not match. So it just prints 1000 / 1000

1 Like
6

Hey thanks for trying to help!
The exact code is beneath the spoiler button

Yes, well, I already know that :slight_smile:

This is correct, but here, and with look rotation in general, 3x3s are guaranteed to be orthogonal, as given by orthogonal basis vectors, which is the special condition where this must work.

That said – and thanks for highlighting the special conditions – I think I’ve found the problem in my analysis. It was the wrong assumption, as always.

I wrongly assumed that taking a cross from two unit vectors must be a unit vector. However, this is the case only when the two vectors are orthogonal to each other (which was also part of my assumption), because the magnitude of cross is the area of the parallelogram.

So when I did

var right = cross(up, forward);

This must be normalized because up and forward can be actually non-orthogonal initially even though they’re of unit length. And then the error propagated from there. I’m still to check out the rest of it, but I think this was the key problem in both of my independent analyses (I had two exactly to root out the common issues but yeah).

Np, thanks, I think this is solved anyway. Assumptions killed the developer.

I will probably come back with LookRotation re-implemented.

@Bunny83
Haha that’s awesome, you somehow averted the error I made (I can see it in line 10). I’ve tried that solution (and I’m still using a cleaned up variant of that for testing). I believe it’s from a decompiled source, looks rather uncouth imo.

I’ve tried what you pasted there, and it does work

However

var rgt = Vector3.Cross(up, fwd);
Debug.Log(rgt.magnitude);

this is 0.333333

So I’m guessing what fixed my problem in your code is on the lines 25 and 26 where these get normalized again.
I’m guessing that’s solved then :slight_smile:

7

On the plus side the source for the method is fully exposed.

https://github.com/Unity-Technologies/UnityCsReference/blob/master/Runtime/Export/Math/Matrix4x4.cs#L370

// Creates a rotation matrix. Note: Assumes unit quaternion
public static Matrix4x4 Rotate(Quaternion q)
{
    // Precalculate coordinate products
    float x = q.x * 2.0F;
    float y = q.y * 2.0F;
    float z = q.z * 2.0F;
    float xx = q.x * x;
    float yy = q.y * y;
    float zz = q.z * z;
    float xy = q.x * y;
    float xz = q.x * z;
    float yz = q.y * z;
    float wx = q.w * x;
    float wy = q.w * y;
    float wz = q.w * z;

    // Calculate 3x3 matrix from orthonormal basis
    Matrix4x4 m;
    m.m00 = 1.0f - (yy + zz); m.m10 = xy + wz; m.m20 = xz - wy; m.m30 = 0.0F;
    m.m01 = xy - wz; m.m11 = 1.0f - (xx + zz); m.m21 = yz + wx; m.m31 = 0.0F;
    m.m02 = xz + wy; m.m12 = yz - wx; m.m22 = 1.0f - (xx + yy); m.m32 = 0.0F;
    m.m03 = 0.0F; m.m13 = 0.0F; m.m23 = 0.0F; m.m33 = 1.0F;
    return m;
}
1 Like
8

I found the full matrix on WP, but this is backwards, quat → 3x3 (4x4 in this case).
And I needed 3x3 → quat which is somewhat harder.

You can’t see the official implementation of LookRotation and I think it’s the same for rotation property on Matrix4x4.

Anyway I found the error, it’s described in the above post. As always it was an SS&S error: something subtle and stupid.

Now all of my variants work :smile: (I’ve accrued a dozen)

Anyway I’m happy to have finally completed the quaternion zoo.
Next project: my own 3D engine, but with blackjack and hookers.

3 Likes
9

Here’s the final implementation

/// <summary> Equivalent to <see cref="UnityEngine.Quaternion.LookRotation"/>. 

/// Assumes unit vectors (directions) on input. </summary>
static public Quaternion LookRotation(Vector3 forward, Vector3 up) {
  if(forward == Vector3.zero) return Quaternion.identity;
  if(up == Vector3.zero) up = Vector3.up;
  var right = Vector3.Cross(up, forward).normalized;
  if(right == Vector3.zero) return Quaternion.FromToRotation(Vector3.forward, forward);
  up = Vector3.Cross(forward, right);
  return ortho3x3ToQuat(ref right, ref up, ref forward);
}

[MethodImpl(MethodImplOptions.AggressiveInlining)]
static Quaternion ortho3x3ToQuat(ref Vector3 c0, ref Vector3 c1, ref Vector3 c2) {
  if(c0.x + c1.y + c2.z > 0f)      return quat(3, 1f + c0.x + c1.y + c2.z, c1.z - c2.y, c2.x - c0.z, c0.y - c1.x);
  if(c0.x >= c1.y && c0.x >= c2.z) return quat(0, 1f + c0.x - c1.y - c2.z, c0.y + c1.x, c0.z + c2.x, c1.z - c2.y);
  if(c1.y > c2.z)                  return quat(1, 1f - c0.x + c1.y - c2.z, c1.x + c0.y, c2.y + c1.z, c2.x - c0.z);
                                   return quat(2, 1f - c0.x - c1.y + c2.z, c2.x + c0.z, c2.y + c1.z, c0.y - c1.x);

  [MethodImpl(MethodImplOptions.AggressiveInlining)]
  static Quaternion quat(int index, float t, float a, float b, float c) {
    var s = .5f / MathF.Sqrt(t);
    return index switch {
      3 => new(a * s, b * s, c * s, t * s),
      2 => new(a * s, b * s, t * s, c * s),
      1 => new(a * s, t * s, b * s, c * s),
      _ => new(t * s, a * s, b * s, c * s)
    };
  }
}

Of course it looks funky and unusual, what did you expect from orion? :slight_smile:

It’s quite serviceable. For me these fundamentals are perfect if they’re well-structured and readable, and I will compactify them until a decent ratio of information entropy is achieved. Apart from switch, there is no repetition and almost no boilerplate. I mean just compare it with the standard decompiled snippet:
Monstrosity

private static Quaternion QuaternionLookRotation(Vector3 forward, Vector3 up)
{
  forward.Normalize();

  Vector3 vector = Vector3.Normalize(forward);
  Vector3 vector2 = Vector3.Normalize(Vector3.Cross(up, vector));
  Vector3 vector3 = Vector3.Cross(vector, vector2);
  var m00 = vector2.x;
  var m01 = vector2.y;
  var m02 = vector2.z;
  var m10 = vector3.x;
  var m11 = vector3.y;
  var m12 = vector3.z;
  var m20 = vector.x;
  var m21 = vector.y;
  var m22 = vector.z;


  float num8 = (m00 + m11) + m22;
  var quaternion = new Quaternion();
  if (num8 > 0f)
  {
    var num = (float)System.Math.Sqrt(num8 + 1f);
    quaternion.w = num * 0.5f;
    num = 0.5f / num;
    quaternion.x = (m12 - m21) * num;
    quaternion.y = (m20 - m02) * num;
    quaternion.z = (m01 - m10) * num;
    return quaternion;
  }
  if ((m00 >= m11) && (m00 >= m22))
  {
    var num7 = (float)System.Math.Sqrt(((1f + m00) - m11) - m22);
    var num4 = 0.5f / num7;
    quaternion.x = 0.5f * num7;
    quaternion.y = (m01 + m10) * num4;
    quaternion.z = (m02 + m20) * num4;
    quaternion.w = (m12 - m21) * num4;
    return quaternion;
  }
  if (m11 > m22)
  {
    var num6 = (float)System.Math.Sqrt(((1f + m11) - m00) - m22);
    var num3 = 0.5f / num6;
    quaternion.x = (m10 + m01) * num3;
    quaternion.y = 0.5f * num6;
    quaternion.z = (m21 + m12) * num3;
    quaternion.w = (m20 - m02) * num3;
    return quaternion;
  }
  var num5 = (float)System.Math.Sqrt(((1f + m22) - m00) - m11);
  var num2 = 0.5f / num5;
  quaternion.x = (m20 + m02) * num2;
  quaternion.y = (m21 + m12) * num2;
  quaternion.z = 0.5f * num5;
  quaternion.w = (m01 - m10) * num2;
  return quaternion;
}

I have manually verified each parameter in every line, a couple of times. It needs more testing, but so far it passes:

  • the visual test (I made a live adjustable visual comparison in the editor),
  • cardinal direction tests + edge cases with zero vectors and/or equal vectors, as well as
  • Bunny’s shotgun test cranked up to 10M iterations (numerous times), showing no discrepancies ( @Bunny83 thanks for the code btw, I was lazy to do it, and it helped me a lot).

I wish I knew how to do a better branching strategy. I did another approach that was also correct, but wasn’t outputting exactly the same results numerically (if all components have their signs flipped, that’s the same quaternion, but Bunny’s test would fail). This delineation turned out to be very complicated, and I left it as it was.

Now, performance. This thing is about 1.55x slower than LookRotation. A single call takes around 202 nanoseconds on my machine on average (in Bunny’s test with random targets) as opposed to 130 nanoseconds for Unity’s LookRotation. To be honest, that’s nothing. It’s perfectly usable. Aggressive inlining matters btw.

Also I made sure to decouple 3x3 → quat conversion from the LookRotation business, so it’s repurposable as well. Hopefully the title will make this searchable.

–
Edit: added xmldoc and static modifiers

2 Likes
10

Well, I haven’t compared the result that you get from Unity’s mathematics package. It uses some fancy bitmasking of the signs to construct several masks. It’s really hard to follow. Though almost all of the float to int conversion are done to directly extract the sign bit and use some xor logic on them, nothing more. The branching seems to be done through the masking and swizzling. Have you actually compared the result of the mathematics implementation with the others? I guess it should also work, no? Since the mathematics package uses all those float3, int3 and uint3 types, some operations look simple but in fact involve several operations. Though I’m almost certain that those can be calculated with the SSE or equivalent SIMD vector opcodes depending on the cpu, especially with burst. Though I haven’t looked into that at all.

1 Like
11

Yep, I checked it out, it’s a great repository, but this particular thing is really hard to follow. Even though I have a good high-level understanding of what’s going on, I’ve decided against trying to parse/transcribe it, because I would’ve wasted days and I still couldn’t tell if I made a typo somewhere. It turned out I was right, plowing through the “internet wisdom” was quicker (if somewhat frustrating).

What I meant by ‘branching strategy’ is just a reorientation of the ‘switching graph’ so to say, which encompasses a smarter mathematics on how to select the dominant component. Specifically, this paper explains it in more detail. Basically someone was very angry with the most commonly found implementation out there (that Unity apparently is using as well).

Here’s the code as presented in this paper.

if (m22 < 0)
 {
   if (m00 > m11)
   {
     t = 1 + m00 - m11 - m22;
     q = quat( t, m01+m10, m20+m02, m12-m21 );
   }
   else
   {
     t = 1 - m00 + m11 - m22;
     q = quat( m01+m10, t, m12+m21, m20-m02 );
   }
 }
 else
 {
   if (m00 < -m11)
   {
     t = 1 - m00 - m11 + m22;
     q = quat( m20+m02, m12+m21, t, m01-m10 );
   }
   else
   {
     t = 1 + m00 + m11 + m22;
     q = quat( m12-m21, m20-m02, m01-m10, t );
   }
 }
 q *= 0.5 / Sqrt(t);

It’s more balanced and vastly superior (in my opinion) to the other method, however it will produce ‘flipped’ (edit: or ‘inverted’) quaternions (xyzw * -1) in ways that make this hard to nail down, at least for my limited understanding of the complex algebra. It delimits this abstract space into some “zones,” and these “zones” overlap with what the other method does in unforeseeable ways (as you can imagine these things don’t correlate with Euclidean space at all), even though the final results are consistently mathematically correct (because (x,y,z,w) == (-x,-y,-z,-w) ).

I’m aware this explanation might be complicated, so let’s try this analogy. Suppose that one solution cuts the “spherical problem domain” (let’s pretend) into two at the equatorial plane, then splits the hemispheres into polar quadrants (half-wedges). Well, the other solution takes the top and down equilateral cones (sectors), then splits the sectors in two, and then splits the remaining equatorial belt in four parts. (Of course I don’t think that’s how quaternions work, it’s just a silly analogy that serves to illustrate two different strategies one can pursue to butcher some sphere into eight pieces.)

The point is, in the end you cannot simply flip some signs in one branch, it’s always messed up somewhere (they’re, erm, let’s say anisomorphic??). But I might try to experiment with this more. If I manage to do this, then this method is a better formulation than what Unity offers (and perhaps even similar performance-wise, which is silly given that it’s on the C++ side, I’m guessing the external call itself introduces a significant overhead; edit: or that the greatest cost of the function lies in the C# part which normalizes vectors).

Edit: changed my wording slightly

12

While I’m at it, here’s a quaternion angle axis decomposer.

/// <summary> Numerically stable conversion from a unit quaternion to angle axis representation (in radians). 

/// Identity quaternions are undefined in this context so the function returns <paramref name="defaultAxis"/>
/// instead. </summary>
/// <param name="axis"> Output. Axis of rotation (unit vector). </param>
/// <param name="rad"> Output. Angle of rotation in the interval -pi..+pi </param>
static public void decompose(Quaternion q, out Vector3 axis, out float rad, Vector3 defaultAxis = default) {
  const float EPS = 1E-8f;
  var sth = MathF.Sqrt(q.x * q.x + q.y * q.y + q.z * q.z);
  if(sth < EPS) { rad = 0f; axis = defaultAxis; return; }
  rad = 2f * (q.w < 0f? MathF.Atan2(-sth, -q.w) : MathF.Atan2(sth, q.w));
  axis = new(q.x / sth, q.y / sth, q.z / sth);
}

Passes all tests, with a caveat:
Quaternion identity is numerically undefined because it doesn’t actually do anything, so no rotation is encoded.

You can pass in a default axis to be returned in that case (the angle will be zero). If the argument is omitted, the function returns a zero vector. This can also be used as a tolerant identity detector, although quaternion dot product perhaps does a better job in this regard (edit: and doesn’t involve a square root).

The returned angle (in radians) is in the interval -pi…+pi (I think it’s inclusive on both sides, but I haven’t tested this specifically). So if you use 270 degrees to construct a rotation, you’ll get back -90 (edit: in radians). A quaternion of (1, 0, 0, 0) i.e. will return pi over (1, 0, 0) axis. Likewise (-1, 0, 0, 0) will return pi over (-1, 0, 0).

Here’s a modified “Bunny’s shotgun” (I think I’ll adopt this as a term lol)
test

void test() {
  const float r2d = Mathf.Rad2Deg;
  const float tau = 2f * MathF.PI;

  int success = 0;
  int C = 1_000_000;

  for(int i = 0; i < C; i++) {
    var axis = Random.onUnitSphere;
    var angle = (tau * Random.value) % tau; // [0..360)

    if(angle < 1E-6f) { --i; continue; } // retries no-op

    var q = Quaternion.AngleAxis(angle * r2d, axis);
    decompose(q, out var retAxis, out var retAngle);

    if(retAngle < 0f) retAngle += tau; // brings it back to [0..360)

    if(MathF.Abs(angle - retAngle) > 1E-6f || retAxis != axis) {
      Debug.Log($"Not equal {(angle * r2d).ToString("f6")} - {(retAngle * r2d).ToString("f6")} ; {axis.ToString("f4")} - {retAxis.ToString("f4")}");
    } else {
      success++;
    }
  }

  Debug.Log($"{success} / {C}");
}

Angle is always correct up to 6 decimal places (half the time up to 7). (edit: this is true unless the angle was very close to multiples of 360 including 0, then the precision might drop slightly and the angle can even turn negative with a small error.)

I haven’t tested this against Quaternion.ToAngleAxis. It should work the same (or similar). I had several version of this lying around for some time, but all of them were acos-based and I was stubborn to find a numerically-stable one. The original source is from this fine repository in C++, which was modified slightly. (More info on why atan2(-sin, -cos) is there).

Again, this is mostly for learning purposes. sth in the above code represents the sine of the angle theta. Whereas q.w is the cosine. This allows to use Atan2, which doesn’t suffer from inaccuracies when the angle is too close to the quadrant boundaries. After that the obfuscated angle axis representation unravels.

–
Edit: removed the need for axis re-normalization (in the sense of redoing sqrt; also redid the tests)