`new Quaternion(Up, Down) * Up == Up`, but it should be Down

[Fruitsalad]

Godot version

4.1

System information

Linux 5.15

Issue description

When using the arc-based quaternion constructor very specifically with the arguments Vector3.Up & Vector3.Down, the resulting quaternion is wrong.

Steps to reproduce

Sorry that it's in C# instead of GDScript, but I don't have much experience with GDScript.

// These assertions are fine:
assert (new Quaternion(Vector3.Left, Vector3.Right) * Vector3.Left == Vector3.Right);
assert (new Quaternion(Vector3.Right, Vector3.Left) * Vector3.Right == Vector3.Left);
assert (new Quaternion(Vector3.Forward, Vector3.Back) * Vector3.Forward == Vector3.Back);
assert (new Quaternion(Vector3.Back, Vector3.Forward) * Vector3.Back == Vector3.Forward);

 // These assertions will fail:
assert (new Quaternion(Vector3.Down, Vector3.Up) * Vector3.Down == Vector3.Up);  // Result is actually Vector3.Down
assert (new Quaternion(Vector3.Up, Vector3.Down) * Vector3.Up == Vector3.Down);  // Result is actually Vector3.Up

[AThousandShips]

What are you basing what it should be on? The other results flipping the vector?

[Fruitsalad]

Based on the description in the documentation, I assumed when an arc on a sphere goes from one side to the other, it should be a 180° rotation. Moreover new Quaternion(new Vector3(0.01f, 0.99f, 0).Normalized(), Vector3.Down).Normalized() * Vector3.Up, gives (0.010114331, -0.99994874, 0) (so almost Down), which seems consistent with what I'd expect it to do. I will admit, however, I'm not very knowledgeable about the workings of quaternions and I might simply be misunderstanding what the function is supposed to do.

[AThousandShips]

For any case where the two vectors are opposite it gives the quaternion (0, 1, 0, 0), which flips through the Y-axis, the case with the vectors being up/down it might break down, but unsure how to handle that seemingly degenerate case

[Fruitsalad]

Again I don't know much about quaternions, but I'd guess flipping either the x-axis or the z-axis would work fine, as long as it's a 180° rotation.

[AThousandShips]

Yes but how do you pick that? I'm not sure this particular use of this constructor is reliable, but I'm not fully familiar with the system

The method for constructing it breaks down when the vectors are parallel, which requires to just pick some value instead

[Fruitsalad]

As far as I can tell there's infinitely many possibilities for the shortest arc between the top and the bottom, so any arbitrary result is fine as long as it's one of the valid results.

[AThousandShips]

But again, how to decide when this is the case? I think this should be documented as I suspect it's quite complex to try to add tests for this specific case

[Fruitsalad]

I believe if the exterior product (a.X * b.Y - a.Y * b.X) of two vectors is 0, they are collinear (by which I mean: facing the same direction or the opposite directions — but I'm not a mathematician so I might be using the wrong word here). I'm not sure how to decide on a result in the general case though.

[AThousandShips]

We check it by seeing if the two vectors are parallel and opposite, tested with the scalar product

My question is: Is this case degenerate? I'd suspect it is and it's not well defined, and the question is if it is worth it to ensure it gives one of the "valid" cases or just return a degenerate case (being clarified in the documentation of course)

Note that it also breaks down in a number of other cases when the input vectors are not along the axes, so that's much more complicated

For example, how would you pick the result for the input vector: Vector3(1, 1, 1).normalized() and the opposite of the same?

[Fruitsalad]

Asking for a 180° rotation seems fairly common for me, so regardless of whether or not it's technically degenerate (which I can't answer), it would seem desirable to me to produce some result that meets the programmer's expectations. Although I can understand that anything to do with quaternions is very difficult and I don't mind this going unresolved for some time. I've already implemented a work-around for my own project.

Just to give a use case, for me it was the following code that triggered this edge case:

void update_transform() {
    var rotation = new Quaternion(Vector3.Up, _normal);
    tf = Rotation(rotation) * Translation(new(0, distance, 0));
}

Which worked quite well except when _normal was (0,-1,0). (With Rotation and Translation producing transformation matrices)

[AThousandShips]

180 degree rotation around what axis though?

[Fruitsalad]

I think any axis should do, as long as it's the right kind of 180° rotation. Again there's infinitely many right answers in this situation, and as long as it's any of those, it seems to me that it should be okay.

[AThousandShips]

But how do you find it?

[Fruitsalad]

It seems to me that any axis that is orthogonal to the two given parallel vectors would be okay, so you could use a cross product to find it.

[AThousandShips]

That is true, the question is if it is worth the several added operations to find this since it is still undefined, there's no such thing as the "shortest arc" in this case

[Fruitsalad]

I thought you said the engine already checks for the case in which the vectors are parallel, in which case it seems to me this isn't really adding any further overhead except in the edge case where the vectors are indeed parallel — and that's the specific case that (in my opinion) currently isn't being handled particularly well. I don't really see any significant performance downside.

[AThousandShips]

Well you'll have two, possibly three, cross products, and an normalization, which is quite a lot more

In total:

• 13 multiplications
• 3 comparisons (for checking if the first cross product is zero)
• One square root

[Fruitsalad]

But it seems to me the added operations would specifically only be in the code branch of the parallel edge case (which doesn't work all that well), and wouldn't affect the normal cases. I might be misunderstanding, though.

[AThousandShips]

It won't, but that doesn't remove the performance increase here

I think there needs to be some argument in favor of this case being a case that isn't degenerate and is worth actually finding a universal fix for

It's really in my opinion a question of finding a less wrong answer

[Fruitsalad]

I think taking a little extra time to produce a good result beats getting a wrong result very quickly.

Edit: I'll await further opinions, nevertheless thank you for your time, especially since it's a volunteer position.

[logzero]

Check if input vectors a, b are antiparallel, which is trivial if you check for parallel already.

If true create a vector which is guaranteed not to be parallel: c = (abs(a.z)+1, a.x, a.y) would work for example.

Compute rotation axis: r = normalize(cross(a, c))

Create 180 deg quaternion: q = (0, r.x, r.y, r.z)

[kleonc]

For me this is clearly a bug. The from, to arguments being the opposite vectors case is already being properly detected[^1], but the always returned 180 rotation around the Y-axis for such cases is correct only for cases when from.y is zero (to.y is also zero then).

If this constructor wouldn't be supposed to handle the opposite vectors then it shouldn't handle them in any way whatsoever. In such case documenting such unsupported behavior and erroring in every such case could be a solution. But as mentioned, this constructor already handles such cases, just most of them incorrectly.

Hence I think for the opposite vectors the proper thing to do is to make this constructor return a rotation by 180 around any arbitrarily choosen axis orthogonal to the input vectors. And document it clearly; users can always manually handle such case if some specific behavior is desired. Likely makes sense to prefer returning rotation around Y axis (when it's correct), this would not change the behavior for the cases already working correctly.

Also looking at the code it seems like it assumes the input vectors are normalized (then dot product is equivalent to the cosine of the angle between the vectors), but it's not ensured/checked anyhow. This should be changed too. https://github.com/godotengine/godot/blob/16a93563bfd3b02ca0a8f6df2026f3a3217f5571/core/math/quaternion.h#L146-L148

[^1]: Given the input vectors are indeed normalized.

[voidexp]

Documentation does not say anything about the case of opposite vectors returning unexpected results, and since a valid 180-degree rotation is still possible around an arbitrarily chosen axis, why don't just return that? With a note, that slerp-ing such rotation might yield unexpected results, since the axis is chosen at random.

[Fruitsalad]

That was more or less the consensus so far (minus the note kleonc made about preferring rotation around the Y axis). It's really just that no one has implemented a fix yet.

[krisutofu]

Check if input vectors a, b are antiparallel, which is trivial if you check for parallel already.

If true create a vector which is guaranteed not to be parallel: c = (abs(a.z)+1, a.x, a.y) would work for example.

Compute rotation axis: r = normalize(cross(a, c))

Create 180 deg quaternion: q = (0, r.x, r.y, r.z)

Crazy, it can be that simple right? Just use an abs().

When I discovered this issue, I only could come up with a general mathematical formula for perpendicular vectors (hint, it's an anti-symmetric matrix with zero trail) but I tried for many hours to find a simple function which would always produce a perpendicular non-zero vector from any non-zero vector. I think, it's impossible with linear or multilinear maps and a kind of case distinction (such as abs()) is necessary.

The next night before sleep, suddenly the idea plopped up, that the minimum absolute value of an input vector needs to be mapped to the maximum absolute value in the perpendicular vector. My own idea therefore is to cross-produce with an approximate perpendicular vector that is 1 at the minimum absolute value and zero elsewhere, and then normalize. Technically, the 1 doesn't need to correspond to the minimum component, just to one component that is small enough. This allows for short-circuit evaluation.

This idea results in the following replacement code which implements the previous paragraph (sorry, I am too busy to deal with forks and pull requests right now). The arguments are

    bool set_perpendicular_to(const real_t& in_x, const real_t& in_y, const real_t& in_z, real_t& out_x, real_t& out_y, real_t& out_z) noexcept {
            constexpr real_t threshold = 1f / Math::sqrt(3f);
            real_t abs_x = Math::abs(in_x);
            if (abs_x > threshold + (real_t)CMP_EPSILON) return false;

            real_t length = Math::sqrt(1f - abs_x*abs_x);
            out_x = 0;
            out_y = -in_z / length;
            out_z = in_y / length;
            return true;
    }

    // v0 and v1 must be unit vectors
    Quaternion(const Vector3 &v0, const Vector3 &v1) { // Shortest arc.
        Vector3 c = v0.cross(v1);
        real_t d = v0.dot(v1);

        if (d < -1.0f + (real_t)CMP_EPSILON) {
            set_perpendicular_to(v0.x, v0.y, v0.z, x, y, z)
                    || set_perpendicular_to(v0.y, v0.z, v0.x, y, z, x)
                    || set_perpendicular_to(v0.z, v0.x, v0.y, z, x, y);
                        w = 0;
        } else {
            real_t s = Math::sqrt((1.0f + d) * 2.0f);
            real_t rs = 1.0f / s;

            x = c.x * rs;
            y = c.y * rs;
            z = c.z * rs;
            w = s * 0.5f;
        }
    }

The helper function better be private. Even better would be throwing a warning [in the Quaternion constructor], if the arguments are not unit vectors. This implementation requires up to 3 additional comparisons and 3 additional calls to Math::abs() when compared to the standard case but therefore, it uses two multiplications less. I did not compile it but theoretically, it should work.

[Fruitsalad]

@krisutofu I took some time this afternoon to implement your code. I should note that my understanding of the inner workings of quaternions is pretty limited but it doesn't seem like anyone else was going to do it. In general it went well except for one assertion I added to the test cases:

Vector3 left(-1, 0, 0);
Vector3 right(1, 0, 0);
Quaternion left_to_right(left, right);

// When the two vectors lie on the XZ plane, the rotation axis should be the
// Y-axis for backwards compatibility.
Vector3 left_to_right_axis = left_to_right.get_axis();
CHECK(left_to_right_axis.abs().is_equal_approx(up));

As kleonc mentioned, when the parameter vectors are on the XZ plane, the resulting quaternion should rotate on the XZ plane because that's consistent with what the old code does. I briefly tried messing around with the code but I figured it's probably better to just ask. Could you let me know how I could adjust the code so that it stays backward compatible?

Currently left_to_right rotates on the XY plane (so around the Z-axis).

@AThousandShips It's very late but I was reading through the thread this morning and I'd like to apologize for the phrasing "since it's a volunteer position and all" because it sounded pretty dismissive when I read it back. I'm not great at social things and my intuition for which tone a text I'm writing has is occasionally pretty bad, but I had always meant this in the sense of "especially since it's a volunteer position." I really didn't mean to be ungrateful. Sorry for that!

[AThousandShips]

No offence taken, thank you for taking the time to clarify!

[krisutofu]

@krisutofu I took some time this afternoon to implement your code. I should note that my understanding of the inner workings of quaternions is pretty limited but it doesn't seem like anyone else was going to do it. In general it went well except for one assertion I added to the test cases:

Vector3 left(-1, 0, 0);
Vector3 right(1, 0, 0);
Quaternion left_to_right(left, right);

// When the two vectors lie on the XZ plane, the rotation axis should be the
// Y-axis for backwards compatibility.
Vector3 left_to_right_axis = left_to_right.get_axis();
CHECK(left_to_right_axis.abs().is_equal_approx(up));

As kleonc mentioned, when the parameter vectors are on the XZ plane, the resulting quaternion should rotate on the XZ plane because that's consistent with what the old code does. I briefly tried messing around with the code but I figured it's probably better to just ask. Could you let me know how I could adjust the code so that it stays backward compatible?

Currently left_to_right rotates on the XY plane (so around the Z-axis).

@AThousandShips It's very late but I was reading through the thread this morning and I'd like to apologize for the phrasing "since it's a volunteer position and all" because it sounded pretty dismissive when I read it back. I'm not great at social things and my intuition for which tone a text I'm writing has is occasionally pretty bad, but I had always meant this in the sense of "especially since it's a volunteer position." I really didn't mean to be ungrateful. Sorry for that!

Great 😊 . Sorry for my response latency. It went into email Spam. The change you need is quite simple, it only needs to swap two things. I assumed X -> Y -> Z order but you'd like to have X -> Z -> Y order.

Therefore the change is to swap these two lines:

        || set_perpendicular_to(v0.y, v0.z, v0.x, y, z, x)
        || set_perpendicular_to(v0.z, v0.x, v0.y, z, x, y);

to

        || set_perpendicular_to(v0.z, v0.x, v0.y, z, x, y)
        || set_perpendicular_to(v0.y, v0.z, v0.x, y, z, x);

XZ has a positive normal that faces downwards (-Y), given that X precedes Z (as it is implemented in the vector comparison). Therefore it gives a -Y as rotation axis (when v0 is Right and v1 is Left). If you still need that to be Y instead of -Y, you can swap the minus sign in the helper function which then gives the negated cross product.

With negated normal, it looks like this:

        out_y = in_z / length;
        out_z = -in_y / length;

Thank you for taking the time for the implementation. We will be grateful.

EDIT: If it relieves you, it's not just you. My communication and social skills are bad in my experience. It takes me a lot of time to find more precise natural language expressions and sometimes I am overthinking (emotionally and such and this makes me hide). As non-native speaker I may notice my phrasing issues less often. But, we are always seeking the benefit of the community, that's the good thing.

[Fruitsalad]

@krisutofu I thought your reply was plenty fast, there's really no need to apologize. I'm not sure if flipping the axis makes a significant difference, both directions seem to work equally well, but based on the old code I guessed that it used to always give a +Y rotation axis, so I've nevertheless added the change for flipping the axis. Thank you for the code and the kind words!