Self Shadow

An online shadow of myself

Perpendicular Possibilities

| Comments

Figure 1: Major axes for original (left), swizzle (mid) and perpendicular (right) vectors

Introduction

Two months ago, there was a question (and subsequent discussion) on Twitter as to how to go about generating a perpendicular unit vector, preferably without branching. It seemed about time that I finally post something more complete on the subject, since there are various ways to go about doing this, as well as a few traps awaiting the unwary programmer.

Solution Quartet

Here are four options with various trade-offs. If you happen to know of any others, by all means let me know and I’ll update this post.

Note: in all of the following approaches, normalisation is left as an optional post-processing step.

Quick ‘n’ Dirty

A quick hack involves taking the cross product of the original unit vector – let’s call it $\mathbf{u}(x, y, z)$ – with a fixed ‘up’ axis, e.g. $(0, 1, 0)$, and then normalising. A problem here is that if the two vectors are very close – or equally, pointing directly away from each other – then the result will be a degenerate vector. However, it’s still a reasonable approach in the context of a camera, if the view direction can be restricted to guard against this. A general solution in this situation is to fall back to an alternative axis:

1
2
3
4
5
6
7
8
9
float3 perp_quick(float3 u)
{
    float3 v;
    if (abs(u.y) < 0.99)          // abs(dot(u, UP)), somewhat arbitrary epsilon
        v = float3(-u.z, 0, u.x); // cross(u, UP)
    else
        v = float3(0, u.z, -u.y); // cross(u, RIGHT)
    return v;
}

Listing 1: A quick way to generate a perpendicular vector

Hughes-Möller

In a neat little Journal of Graphics Tools paper [1], Hughes and Möller proposed a more systematic approach to computing a perpendicular vector. Here’s the heart of it:

Or, as the paper also states: “Take the smallest entry (in absolute value) of $\mathbf{u}$ and set it to zero; swap the other two entries and negate the first of them”.

Figure 2: Distribution of v over the sphere

However, there’s a problem with this as written: it doesn’t handle cases where multiple components are the smallest, such as $(0, 0, 1)$! I hit this a few years back when I needed to generate an orthonormal basis for some offline geometry processing, and it’s easily remedied by replacing $<$ with $\leq$. Here’s a corrected version in code form:

1
2
3
4
5
6
7
8
9
10
11
12
float3 perp_hm(float3 u)
{
    float3 a = abs(u);
    float3 v;
    if (a.x <= a.y && a.x <= a.z)
        v = float3(0, -u.z, u.y);
    else if (a.y <= a.x && a.y <= a.z)
        v = float3(-u.z, 0, u.x);
    else
        v = float3(-u.y, u.x, 0);
    return v;
}

Listing 2: Hughes-Möller perpendicular vector generation

Stark

More recently, Michael M. Stark suggested some improvements to the Hughes-Möller approach [2]. Firstly, his choice of permuted vectors is almost the same – differing only in signs – but even easier to remember:

In plain English: the perpendicular vector $\mathbf{\bar{v}}$ is found by taking the cross product of $\mathbf{u}$ with the axis of its smallest component. (Note: the same care is needed when multiple components are the smallest.)

Figure 3 visualises this intermediate ‘swizzle’ vector over the sphere:

Figure 3: Intermediate swizzle vector

Secondly, Michael also provides a branch-free implementation:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
float3 perp_stark(float3 u)
{
    float3 a = abs(u);
    uint uyx = SIGNBIT(a.x - a.y);
    uint uzx = SIGNBIT(a.x - a.z);
    uint uzy = SIGNBIT(a.y - a.z);

    uint xm = uyx & uzx;
    uint ym = (1^xm) & uzy;
    uint zm = 1^(xm & ym);

    float3 v = cross(u, float3(xm, ym, zm));
    return v;
}

Listing 3: Branch-free perpendicular vector generation

I know what you’re thinking, there’s a problem here too: it should be zm = 1^(xm | ym)! Although still robust, the effect of this error is that the nice property of even, symmetrical distribution over the sphere is lost:

Figure 4: Broken symmetry in the perpendicular vector
(that’s 7 years bad luck!)

XNAMath

Finally, another branch-free solution is provided in the form of XMVector3Orthogonal, which is part of the XNAMath library. Here’s the actual code taken from the DirectX SDK (June 2010):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
XMFINLINE XMVECTOR XMVector3Orthogonal(FXMVECTOR V)
{
    XMVECTOR NegativeV;
    XMVECTOR Z, YZYY;
    XMVECTOR ZIsNegative, YZYYIsNegative;
    XMVECTOR S, D;
    XMVECTOR R0, R1;
    XMVECTOR Select;
    XMVECTOR Zero;
    XMVECTOR Result;
    static CONST XMVECTORU32 Permute1X0X0X0X =
        {XM_PERMUTE_1X, XM_PERMUTE_0X, XM_PERMUTE_0X, XM_PERMUTE_0X};
    static CONST XMVECTORU32 Permute0Y0Z0Y0Y =
        {XM_PERMUTE_0Y, XM_PERMUTE_0Z, XM_PERMUTE_0Y, XM_PERMUTE_0Y};

    Zero = XMVectorZero();
    Z = XMVectorSplatZ(V);
    YZYY = XMVectorPermute(V, V, Permute0Y0Z0Y0Y.v);

    NegativeV = XMVectorSubtract(Zero, V);

    ZIsNegative = XMVectorLess(Z, Zero);
    YZYYIsNegative = XMVectorLess(YZYY, Zero);

    S = XMVectorAdd(YZYY, Z);
    D = XMVectorSubtract(YZYY, Z);

    Select = XMVectorEqualInt(ZIsNegative, YZYYIsNegative);

    R0 = XMVectorPermute(NegativeV, S, Permute1X0X0X0X.v);
    R1 = XMVectorPermute(V, D, Permute1X0X0X0X.v);

    Result = XMVectorSelect(R1, R0, Select);

    return Result;
}

Listing 4: XNAMath’s semi-vectorised method

Let me save you the trouble of parsing this fully (or consider it an exercise for later); if you boil it down, what you’re effectively left with is:

I’ve failed, thus far, to pinpoint the origin of or thought process behind this approach. That said, some insight can be gained from visualising the resulting vectors:

Figure 5: Covering one’s axis

Their maximum component is in the x axis, except close to the +/-ve x poles. Essentially, Microsoft’s solution ensures robustness without concern for distribution, much like the initial ‘quick’ approach.

Performance

I haven’t benchmarked these implementations, since in cases where I’ve needed to generate perpendicular vectors, absolute speed wasn’t important or the call frequency was vanishingly small. Even in performance-critical situations, it really depends on what properties/restrictions you can live with and your target architecture(s). Still, I can’t help but think that XMVector3Orthogonal is doing a little bit more than it needs to, so maybe there’s cause to revisit this subject at a later date.

Conclusion

I hope you’ve learnt something about generating perpendicular vectors, or that I’ve at least made you aware of some of the minor issues in previous work on the subject. On that note, if you spot any new errors here, please let me know!

References

[1] Hughes, J. F., Möller, T., “Building an Orthonormal Basis from a Unit Vector”, Journal of Graphics Tools 4:4 (1999), 33-35.
[2] Stark, M. M., “Efficient Construction of Perpendicular Vectors without Branching”, Journal of Graphics Tools 14:1 (2009), 55-61.

Comments