Quaternions: Handle 180 degree arcs in the shortest arc constructor better

[Fruitsalad]

This PR resolves #80249, where the problem was that the shortest arc quaternion constructor didn't work well when asked to make a 180° quaternion (which is notably ambiguous because there's no single shortest arc for a 180° rotation). Previously only the case where both input vectors lay on the XZ plane was handled correctly, and any other pair of input vectors would give a completely incorrect result. Now any pair of vectors should give a sensible result.

The main part of the implementation is by @krisutofu, but he did not have the time to do the pull request. I've also added a unit test, although it's fairly basic due to my limited understanding of the inner workings of quaternions.

Some notes:

• The implementation required a helper function, but I wasn't sure whether it would be better to make it a C-style static function or a C++-style private function (which are both kind of the same concept). I used a static function but let me know if a private function is preferred.
• The big try_set_perpendicular_to line is so long it's kind of ugly, but that's clang-tidy's handywork.

[Fruitsalad]

I've adjusted the code based on your feedback just now.

[Mickeon]

I wish this made some sense to me but it's fairly difficult to unravel.

What I can see for sure is that this PR makes the Quaternion constructor a whole lot more complicated, and given how significant Quaternions are for the engine, it's worth a pretty close look.

[AThousandShips]

Indeed it does, and IMO the question is if it's worth it for a degenerate case, there isn't any one valid result of this case, so I'm not convinced spending more performance to generate a different, not unique result to a question that doesn't have an answer

It's asking "which is the shortest path from the top of a sphere to the bottom of a sphere"

My opinion is that this should be documented instead

[Fruitsalad]

I've marked it unlikely. In my opinion you can pose the question as "what is a shortest arc?". Just like in a shortest path algorithm, there doesn't necessarily need to be a single correct answer, and in case of multiple correct answers just giving one of them is OK.

[AThousandShips]

Yes but all of them are equally wrong as well as equally right, so making it a different, not right not wrong answer needs to be justified if it's a more costly one, I don't know that it's reasonable to rely on the result of this operation, it doesn't give a useful answer in general, so IMO no one should be using it or need it to be differently right/wrong

The only difference between this particular case and the other ones (unmodified) is that this one doesn't process two specific vectors the same way, but all the other vectors are still essentially randomly transformed, why is that specific part relevant?

I.e.: When would or even should anyone rely on the result of this particular configuration? Is it safe or reasonable to use it? Is it even useful or meaningful? What actual relevant uses would there be for this case?

Edit (to not be too noisy): If you want to turn your spaceship from pointing slightly above the horizon, while also being slightly turned to the side, to facing the other way, will this construction retain reliably the up and side directions? This nominally makes the forward direction flipped, but every other direction should still be just unreliable. But I'm not willing to die on this hill, but given what I know about quaternions I'm very unconvinced that the result of this, "fixed" or otherwise, yields results anyone should rely on in the first place

[Fruitsalad]

I do agree it's worth documenting, since there's no single shortest arc here, but I also think this PR is frankly an upgrade over the old code. I don't really agree with all of the arcs being "equally wrong" since that would imply it's not a valid shortest arc, the way I see it they're just all "equally right". For my argument I'll make a plea to popularity though: About six/seven people ran into the same problem and agreed with the original issue in the past ten months, so it does seem to me that people expect it to work this way.

I can't really give good use cases right now since I've only needed this once, which was about a year ago. Back then I was generating a quaternion to interpolate between two direction vectors, I believe, and of course the old/current way in which 180° rotation is handled in the arc constructor lead to an unexpected bug. The new way is arguably not-entirely-correct but I think the old way is blatantly incorrect. I will note that either documenting it or adding an error message for 180° rotations also works and wouldn't have lead to a problem in my case. Maybe someone else can chime in on use cases.

Edit: I don't actually know about the up vectors question. I think if you multiply the arc-based quaternion by the spaceship's old quaternion it might preserve a sensible up vector? I'm not confident about that, though.

Edit2: Actually you're probably right about the up vectors matter. It might just as well end up upside-down as it might end up right side up in the 180° case.

[AThousandShips]

I'd say that if it wasn't for compatibility we should just error out and return an empty quat in this case, as I'd say the return value currently (unmodified) is just a "well, I have to answer something", just like with looking_at which errors in degenerate cases

But it's unfortunately not really a viable option here I'd say, though I'm still not convinced the opposite case is something is or should be using so arguably changing it to error might be justified

[Mickeon]

It's worth noting that even with look_at/looking_at there have been "requests" to figure out a best guess, because the current behavior is seen as annoying, even if it's mathematically impossible.
This is at least from what I've seen in casual conversation, or for example #57457 and #53793 (comment)

[Fruitsalad]

I think keeping the old code but adding a warning to the 180° branch is arguably backwards compatible, since it's still allowed just "strongly suggested not to do that". Of course that would be obnoxious for someone happily using the 180° branch that there isn't much difference between a warning and an error in practice, but the current 180° code really only works in the specific case of rotation around the Y axis, in which case you might as well inform people to use an angle-axis constructor.

[AThousandShips]

A warning might be a bit noisy, but with a warn once maybe it would be more useful, and adding details to the documentation

[krisutofu]

I wish this made some sense to me but it's fairly difficult to unravel.

What I can see for sure is that this PR makes the Quaternion constructor a whole lot more complicated, and given how significant Quaternions are for the engine, it's worth a pretty close look.

I do not agree on the argument of being complicated. It's the minimal useful solution that is in theory about as performant as the standard case.

[deleted lengthy irrelevant text]

I read from your posts, that you actually want an explanation. Here it is.

Explanation

It selects a bitangent (or tangent) vector in the rotation plane to compute the rotation axis between two opposite vectors (i.e. quaternion angle = 90°). The bitangent vector is selected from (1, 0, 0), (0, 1, 0) or (0, 0, 1) because these produce the computationally simplest rotation axes.

Perpendicular vectors (with small scalar product) have the property that they have a maximum value in the component where the other vector has a minimum value. Therefore the code tests whether the absolute value of the tested input component is too large for a minimum (or small) component and tries the next axis if the answer is yes. (It excludes the maximum vector component without visiting all components!)

The threshold explanation is, Vector3(1,1,1)/sqrt(3) is the unit vector with 3 equal components. (Which must satisfy the first test.) Therefore, it is the mathematical upper bound of the minimum absolute value of unit vector components.

If an axis is selected, the bitangent vector is implicitly cross produced with the input vector (tangent), and normalized, to obtain the final rotation axis.

The main problem is, there is no linear mapping of unit vectors to perpendicular vectors. There are maybe continous mappings but I assume they have more computational effort. (There is no continous mapping.)

In the function call _try_set_perpendicular_to(…), the rotation of the vector components x, y and z allows the selection of the bitangent and the rotation axis to be applied to a different axis of the original vector. It avoids the bad idea of duplicating code into nested ifs.

Philosphy

Good software doesn't focus on code size or the simplicity of the implementation but usability and utility for the user because that is the most important criterion for choosing software. If you don't want to remove this quaternion constructor entirely (making breaking changes), you need any useful solution to handle these edge cases.

[krisutofu]

I'd be asking for only three possible changes in the PR.

1. Performance improvements:

• Is the 1 / sqrt(3) + CMP_EPSILON treated as constant at runtime? If not it would lead to a relevant loss of Speedup when it is evaluated at runtime. We could make use of a static variable (e.g. with a name like maximum_component_threshold) or it could be replaced by the constant value as a float literal (with a comment of the formula).

  static real_t _maximum_component_threshold = 1f / Math::sqrt(3f) + (real_t) CMP_EPSILON;
  …
      if (abs_x > Quaternion::_maximum_component_threshold) return false;

• I figured out, it's faster to compute the inverse square root directly (with a compiler intrinsic or buildin) into one variable and multiplying it, instead of computing the square root and dividing it two times. If you want best performance (instead of highly accurate rotations), it could attempt to use compiler intrinsics/builtins which may compile to an inverse-square-root CPU instruction for certain CPU targets. The result is probably notably worse if the radius becomes large.

2. Better game developer experience: It would be better to add a warning, when users plug in a non-unit vectors. Without unit vectors it will not work as expected, it can easily produce a NaN quaternion.

3) If it doesn't matter whether the right vector uses +Y or -Y as rotation axis, it would be preferable to -Y and use the current tested axis as tangent instead of as bitangent.

r_out_x = 0;
r_out_y = -p_in_z;
r_out_z = p_in_y;

The reason is, +Y axis rotates the right vector counterclockwise in the XZ plane (using the formula v' = p · v · p⁻¹). This matters if users decide to linearly interpolate the quaternion rotation. Clockwise direction is positive in XZ whereas it is negative in XY when we use a right-hand side coordinate system as Godot does in 3D space. The convention would be positive direction for positive angles.

[TokageItLab]

In my opinion, I rather agree that this works. The reason is that Quaternion's 180 degree slerp is valid and often used.

However, I do have a concern about the inconsistency about look_at not working #93653 (comment). This is because it is the same error due to the fact that the choices regarding parallel vectors are infinite, as this PR is trying to solve.

If a method could be defined in the Math class or somewhere else to solve the parallel problem so that it could be shared with other classes and solve the look_at problem as well, I think we would be fine.

Also, there is a slight concern about consistency with slerp. So some additional testing may be needed to ensure that no new issues of inconsistency between slerp and generated arc are raised there.

[krisutofu]

In my opinion, I rather agree that this works. The reason is that Quaternion's 180 degree slerp is valid and often used.

However, I do have a concern about the inconsistency about look_at not working #93653 (comment). This is because it is the same error due to the fact that the choices regarding parallel vectors are infinite, as this PR is trying to solve.

If a method could be defined in the Math class or somewhere else to solve the parallel problem so that it could be shared with other classes and solve the look_at problem as well, I think we would be fine.

Also, there is a slight concern about consistency with slerp. So some additional testing may be needed to ensure that no new issues of inconsistency between slerp and generated arc are raised there.

Oh I see. The problems look similar but they are somewhat different. I've "looked at" the look_at() issue you mentioned.

EDIT: It took me some time to realize the use case. The solution of the look_at problem is very simple, just using a default value for the X axis when the cross product got zero. The default value needs to be obtained from the current transform of the Node3D (EDIT: multiplied with the sign between target and start direction; this way, the basis change remains minimal as desired). But if there is no transform, if no default value is provided (and this is the only implemented case at the moment), it needs to select a vector that is perpendicular to the target direction. Since there is no correct or wrong answer, as long as it is perpendicular, it satisfies the specification.

[Fruitsalad]

@krisutofu I actually forgot about 1.0f / Math::sqrt(3.0f) — it was a constexpr in your original code but sqrt isn't constexpr, so yesterday I fixed it "for the time being" by making it const and then forgot about it entirely. Sorry about that. In any case I'll await whatever conclusion to maintainers come to before making any further changes to the code.

[krisutofu]

Btw, there is an alternative to write it like this:

…
static void ::_set_perpendicular_to(real_t p_x_squared, real_t p_y, real_t p_z, real_t p_square_length const real_t &r_out_x, const real_t &r_out_y, const real_t &r_out_z) {
    auto inverse_length = 1f / Math::sqrt(p_square_length - p_x_squared);
    r_out_x = 0;
    r_out_y = - p_z * inverse_length; 
    r_out_z = p_y * inverse_length;
}

…
if (auto x_squared = p_v.x * p_v.x; x_squared <= (1f/3f) + (real_t)CMP_EPSILON) {
    _set_perpendicular_to(x_squared, p_v.y, p_v.z, 1f, x, y, z);
}
else if (auto y_squared = p_v.y * p_v.y; y_squred < (1f/3f) + (real_t)CMP_EPSILON) {
    _set_perpendicular_to(y_squared, p_v.z, p_v.x, 1f, y, z, x);
} else {
    _set_perpendicular_to(p_v.z * p_v.z, p_v.x, p_v.y, 1f, z, x, y);
}

Maybe, this is easier to read and optimize by the compiler.

[fire]

Here are the concerns:

• Animation team reviewed
• We should fix this. The reason is that Quaternion's 180 degree slerp is valid and often used.
• Reduce code duplication in look_at and other places in Godot Engine
• Check possible inconsistency with slerp and generated arc.
• Non-animation team concerns involve performance which would require before and after benchmarking on hundreds of skeletons or whatever uses this use case.

[TokageItLab]

@fire Since the solver's calculation is only done when the axes are parallel in look_at() or quaternion constructor by arc, whatever the algorithm is, I think that the performance is not a concern in many cases.

But regardless, the algorithm seems to be a bit more optimizable. Since Godot does not optimize for division in the default compilation settings, it is recommended that multiple divisions be replaced with multiplication or cached if precision is not an issue.

[krisutofu]

Hey guys, I just came back thinking about the mapping. Now I understand the criticism and I apologize for my misunderstanding. I begin to see that it is "overcomplicated".

At the moment, it works like this:
Vectors with x between -1/sqrt(3) and 1/sqrt(3) are rotated through the X axis. Else, vectors with z between -1/sqrt(3) and 1/sqrt(3) are rotated through the Z axis. All remaining vectors are rotated through Y.

I admit, the threshold could also be different. We could get rid of the rotations through the Y axis by setting the threshold to 1/sqrt(2) or bigger, so that only X and Z are remaining.

---

Then I noticed, there is a simple almost-continous mapping which is typically used in computer graphics applications. (It makes me feel stupid.) And it is likely the solution that produces most predictable results for users.

We could define one single axis where all vectors should be rotated through, e.g. -Z (the forwards-looking direction). I think, -Z is good because it is useful from a player perspective (rotations will pass through the center of the camera). Then we simply handle the singularity case.

This means, the normalized cross product vector × (0, 0, -1) = (0, 0, 1) × vector is computed like the PR does at the moment. If the cross product is too close to zero, it would use (0, 1, 0) as rotation axis (default rotation on the ground). This would fit the currently implemented behaviour even if there is no continuity in the singularity.

It looks like this could have been the original code intent but somehow, someone forgot to add the cross product into it.

The lack of continuity in one point is still bad, e.g. when the Quaternions are constructed in a sequence, but at least it satisfies the current specification.

The optimal solution lets users optionally pass more information to the constructor to get a situation-specific better result than the default one. If they fail to provide a useful value, it's the fault of the user and they can recover from it by fixing their input value.

For consistency, default behaviour could be put into a function that is reused across the code base. So that users know, everytime, when the rotation gets ambigous and no more information is provided, it tries to rotate through -Z by default, using +Y as rotation axis if also -Z is ambiguous. This is simple to memorize, simpler than a subdivision of the sphere.

New code suggestion (pseudo C++):

Quaternion::_set_default_rotation_axis(const Vector3 &p_v0) {  // or a public static function
    if (Math::is_zero_approx(p_v0.y) && Math::is_zero_approx(p_v0.x)) {
        r_out_x = 0;
        r_out_y = 1;
        r_out_z = 0;
        return;
    }
    real_t inverse_length = 1f / Math::sqrt(1f - p_v0.z * p_v0.z);
    r_out_x = -p_v0.y * inverse_length;
    r_out_y = p_v0.x * inverse_length;
    r_out_z = 0f;
}

// documentation: p_via is a vector in the rotation arc and is only used if p_from and p_target are 180° apart.
Quaternion::Quaternion(const Basis &p_from, const Vector3 &p_target, const Vector3 &p_via) {
…
if (…) {  // edge case, p_from.z.cross(p_target) is zero
    auto rotation_axis = p_from.cross(p_via);
    if (rotation_axis.is_zero_approx())
        this->_set_default_rotation_axis(p_from);   // blame the user, if this is not the desired direction
    else {
        rotation.normalize();
        [r_out_x, r_out_y, r_out_z] = rotation_axis;
    }
    return;
}
…
}

Quaternion::Quaternion(const Vector3 &v0, const Vector3 &v1) {
…
if (…) {  // edge case, assume rotation through -Z with +Y as fallback rotation axis
    this->_set_default_rotation_axis(p_v0);
    return;
}
…
}

It would be interesting, if we could specify a rotation which goes through 3 vectors but this would mean we need an encoding for arcs over 180°. Quaternions only encode the shortest arc.

[krisutofu]

Apologies for writing again.

I found there is a main problem for Godot users, with the idea of rotating through a static fixed axis. It can produce rotations through +Y or -Y axis and the rotation arcs will bend strangely.

User centered rotation design

I wondered what would be the most useful definition for users?

If we want the player, the object, the spaceship, the …, to have a desirable basis after the rotation, we could say, the changes into the Y direction (of any basis vector) should be minimal. This assumes, the players, the objects, the spaceships, are grounded. It is very rare in games to see 180° rotations which leave the object upside down or just lopsided.

The answer here is, to choose a dynamic rotation bitangent. Then compute the rotation axis. We want to prevent changes in Y components of vectors, so this bitangent should always be located in the XZ plane. The minimal extent of Y movement is obtained by making this bitangent the bisectrix of start and target vector (i.e. 90° away from start and end vector).

One unknown is left. We know, there are two bisectrices for an 180° angle. Should it only rotate counter clockwise, clockwise, or should it switch between clockwise and counter clockwise when it transitions between left and right half space or between front and back half space? Either choice of them is fine I believe. It won't affect the basis after the rotation and will not lead to artificial arcs.

Choosing one of two possible rotation directions

One possible choice is "only rotate clockwise" because clockwise rotation in XZ is the positive rotation direction (and we are used to clockwise rotations anyways, it is reminiscent of the flow direction of text). But is this the user-prefered choice?

A criterion for this choice is to minimize the coding effort for Godot users. It needs to cater to a general case or the most common case. Let's look at some cases for "only rotate clockwise".

Case 1:
a game uses only unidirectional rotations, say, only clockwise rotations and never counter-clockwise rotations. Success.

Case 2:
A game uses only counter-clockwise rotations. The user needs to handle the special case of 180° with own code or prevent it from happening. (This issue can be solved by letting users specify an optional "outside" vector, I'll explain below.)

Case 3:
a game uses clockwise directions in the right half space of the vector's coordinate system and counter-clockwise rotations in the left half space. Half the cases are fine, other half of cases are weird. I don't know, if this is more terrible than case 2. This is also solved by specifying an "outside" vector which the user chooses based on the start vector. This is very simple. The case distinction is made by the user, the Quaternion constructor then does the rest of computing the rotation axis.

Case 4:
The user wants a rotation against our expectation of minimal Y rotation or movement but they still want the basis to be plausible after the rotation. Example: flip the vehicle/object into the opposite direction, through the air, such that it lands on the ground again. Actually this means two chained rotations. First, it needs to rotate along the original forwards vector as rotation axis, then rotating the whole result basis with this Quaternion constructor. It's not a single quaternion but a quaternion function over the product of two quaternions which is applied to the basis at animation start.

It can be split into two separate rotations where the 2nd rotatation might use this Quaternion constructor. If the Quaternion constructor is called with Vector3.up to Vector3.down, the right vector could be chosen automatically as default value (as a special case because Vector3.up is unavailable for this case). I would guess, it is more common to rotate around the right vector as axis, then the left vector. When this is implemented, no optional vector needs to passed.

performance

If we assume, that it is the most reasonable solution, that is worthy to provide to users, we need to accept a few additional computational overhead. The usability and utilty will create more value for the user than a little performance improvement.

Algorithm

The outside-vector describes the vector that is always excluded from (outside of) the rotated plane. In other terms, it lies in one plane with start vector, end vector and rotation axis. If this isn't possible, because it is a scaled multiple of the start or end vector, then assume Vector3.up is used as outside vector (the sensible default case). If this is insufficient, use Vector3.right as last resort. (This is useful for case 4, where.)

The final rotation axis of the quaternion is computed with the normalized cross product (v0 × outside).normalized() × v0 (note, the tangent is (v0 × outside).normalized()). It assumes, v0 is the start vector, a unit vector.

Code

Demonstration (bloating performance optimization).

private:
// normalized((x,y,z) × (1, 0, 0)) × (x,y,z)  if (x,y,z) != (1, 0, 0)
static void Quaternion::_set_perpendicular_using_x(real_t p_x, real_t p_y, real_t p_z, real_t &r_x, real_t &r_y, real_t &r_z) {
    real_t inverse_length = 1f / Math::sqrt(1f - p_x * p_x);
    r_x = (p_z * p_z + p_y * p_y) * inverse_length;
    r_y = -p_y * p_x * inverse_length;
    r_z = -p_z * p_x * inverse_length;
}

// normalized((x0, y0, z0) × (0, y1, z1))  if (x0, y0, z0) != (0, y1, z1)
static void Quaternion::_normalized_cross_with_zero_in_x(real_t p_x0, real_t p_y0, real_t p_z0, real_t p_x1, real_t p_y1, real_t p_z1, real_t &r_x, real_t &r_y, real_t &r_z) {
    r_x = (p_y0 * p_z1 - p_z0 * p_y1);
    r_y = -p_x0 * p_z1;
    r_z = p_x0 * p_y1;

     real_t inverse_length = 1f / Math::sqrt(r_x * r_x + r_y * r_y + r_z * r_z);
     r_x *= inverse_length;
     r_y *= inverse_length;
     r_z *= inverse_length;
}

void Quaternion::_set_sign(real_t sign) {
    x *= sign;
    y *= sign;
    z *= sign;
}

// choose the vector, that is perpendicular to p_v0 and in one plane with p_v0 and p_outside
// else try p_outside = (0, 1, 0), else p_outside = (1, 0, 0).
void Quaternion::_set_perpendicular_to(const Vector3 &p_v0, std::optional<Vector3> p_outside) {
    // v0 × outside × v0
    real_t v0_x = p_v0.x;
    real_t v0_y = p_v0.y;
    real_t v0_z = p_v0.z;
    if (p_outside && !p_v0.is_equal_approx(*p_outside)) {
        Vector3 tangent;
        // performance optimization, fastest computation for each case
        if (Math::is_zero_approx(p_outside->y)) {  // y = 1 is the default, therefore a user-specified y = 0 is more likely
            if (Math::is_zero_approx(p_outside->x)) {
                _set_perpendicular_using_x(v0_z, v0_x, v0_y, z, x, y);
                _set_sign(p_outside->z);
                return;
            }
            else if (Math::is_zero_approx(p_outside->z)) {
                _set_perpendicular_using_x(v0_x, v0_y, v0_z, x, y, z);
                _set_sign(p_outside->x);
                return;
            }
            else {
                _normalized_cross_with_zero_in_x(v0_y, v0_z, v0_x, 0f, p_outside->z, p_outside->x, tangent.y, tangent.z, tangent.x);
            }
        } else {
            if (Math::is_zero_approx(p_outside->x)) {
                if (Math::is_zero_approx(p_outside->z)) {
                    _set_perpendicular_using_x(v0_y, v0_z, v0_x, y, v, z);
                    _set_sign(p_outside->y);
                    return;
                } else
                    _normalized_cross_with_zero_in_x(v0_x, v0_y, v0_z, 0f, p_outside->y, p_outside->z, tangent.x, tangent.y, tangent.z);
            } else if (Math::is_zero_approx(p_outside->z)) {
               _normalized_cross_with_zero_in_x(v0_z, v0_x, v0_y, 0f, p_outside->x, p_outside->y, tangent.z, tangent.x, tangent.y);
            } else {
                tangent = p_v0.cross(*p_outside);
                real_t inverse_length = 1f / tangent.length();  // length cannot be zero
                tangent.x *= inverse_length;
                tangent.y *= inverse_length;
                tangent.z *= inverse_length;
            }
        }

        [x, y, z] = tangent.cross(p_v0);
        return;
    }
    // if (0, 1, 0) cannot be used as default for p_outside
    if ( unlikely(Math::is_zero_approx(p_v0.x) && Math::is_zero_approx(p_v0.z)) ) {
        rotate_around_x:
        _set_perpendicular_using_x(v0_x, v0_y, v0_z, x, y, z);
        return;
    }

    rotate_around_y:
    _set_perpendicular_using_x(v0_y, v0_z, v0_x, y, z, x);
}

public:
Quaternion::Quaternion(const Vector3 &p_v0, const Vector3 &p_v1) : Quaternion(p_v0, p_v1, std::nullopt) {}
Quaternion::Quaternion(const Vector3 &p_v0, const Vector3 &p_v1, std::optional<Vector3> p_outside)
{
    …
    if (…) {   // special case
        _set_perpendicular_to(p_v0, p_outside);
        return;
    }
    …
}

If the case-specific optimization is not desired, it can be reduced to the cases with two zeroes in p_outside.

[krisutofu]

I should have known it but I realized, normalize(v0 × outside) × v0 is equivalent to normalized(outside.slide(v0)). This reduces computational effort.

Idea: first try normalize(outside.slide(v0)) else try normalize(vector3(0, 1, 0).slide(v0)) else try normalize(vector3(1, 0, 0).slide(v0)). This computes the rotation axis as unit vector.

The case with two zeros and outside = (1,0,0) stays almost the same. We can just use that y² + z² = 1 - x² so that (y² + z²) / sqrt(1 - x²) = sqrt(1 - x²).

out_x = sqrt(1 - x*x)
temp = 1 / out_x
out_y = -y * x * temp
out_z = -z * x * temp

For the other cases, if outside is set, outside.slide(v0) is computed, checking against zero before normalizing it. If it was zero, it either attempts (0, 1, 0).slide(v0) and normalizes the result if v0 != (0, _, 0) or otherwise it uses (v0.y, 0, 0) . Note that v0 × (1, 0, 0) × v0 = (v0.y, 0, 0) if v0.x == 0 and v0.z == 0.

If no outside is set, it directly continues with (0, 1, 0).