A good way to find a vector perpendicular to another vector?

Question

What is a good way to generate an arbitrary vector that is perpendicular to a vector? Given that the generated vector should be valid (i.e. unless the given vector is (0.0, 0.0, 0.0) - the output for this case should be error)?

I know that there are infinite set of possible vectors (for example see here). What I want is something rigid and bug-free for any input vector.

What I attempted was setting, here p and n are said perpendicular vectors: p[0] = 0.5*(n[0] + n[1]; p[1] = 0.5*(n[0] - n[1]) and then solve for the dot product equation to find p[2] but this is flawed when n[2] is 0. At the moment I can only think of tediously creating all the case check - but there must be better and more elegant solution?

Answer 1

Pick an arbitrary vector that is not collinear with the given vector. For "something rigid", you can have a fixed rule. For instance, pick (1, 0, 0) unless the vector is (x, 0, 0) for some x; otherwise pick (0, 1, 0).¹ Then take the cross product of the arbitrary vector with the input vector. That will yield a perpendicular vector. (This only works in 3D, of course, but that seems to be what you are working in.)

Alternatively (and this works in any dimension), pick the unit vector along the coordinate axis that yields the smallest (in magnitude) dot product with the input vector. Call this unit vector e and the input vector x. Then e − (e • x) x will be perpendicular to x. (It's easy to check that the dot product is zero: assume, without loss of generality, that x is a unit vector. Then (e − (e • x) x) • x = (e • x) − (e • x)(x • x) = (e • x) − (e • x)(1) = 0.)

1_{A better rule would be to pick whichever of (1, 0, 0), (0, 1, 0), or (0, 0, 1) forms the smallest dot product (in magnitude) with the input vector. That will give you better numerical stability.}

Answer 2

Take the cross product of the given vector with a standard basis vector. For best accuracy, take the basis vector which is the less parallel, i.e. which forms the smallest (absolute) dot product. No need of big hammers:

if |Nx| = min(|Nx|, |Ny|, |Nz|) => (0, -Nz, Ny)
if |Ny| = min(|Nx|, |Ny|, |Nz|) => (Nz, 0, -Nx)
if |Nz| = min(|Nx|, |Ny|, |Nz|) => (-Ny, Nx, 0)

Code:

double Ax= abs(N.x), Ay= abs(N.y), Az= abs(N.z); 
if (Ax < Ay)
    P= Ax < Az ? Point(0, -N.z, N.y) : Point(-N.y, N.x, 0);
else
    P= Ay < Az ? Point(N.z, 0, -N.x) : Point(-N.y, N.x, 0);

Answer 3

You could use linear algebra or different trig functions to solve this problem.

Given a vector A(xi,yj,zk) where i,j,k are the unit axis vectors and x,y,z are the magnitudes of that vector you can construct a line by choosing two separate points that exist on the projection of that vector so that points P1(xi1,yj1,zk1) and P2(xi2,yj2,zk2) exist on the line L made by vector A. We know that a line segment of L = P2 - P1 and exists or is parallel to A thus giving you a dot product of 1 and if the dot product is 0 they are Orthogonal or Perpendicular. From these two points you can then find the slope of that line in 3D using the math found in the accepted answer to this post StackExchange:Mathematics and once you have the slope of that line and from basic mathematics; any line that is orthogonal to another their products will = -1. In 2D Euclidean Geometry if a slope is 1/2 its perpendicular will have a slope of -2; so once you find the slope of that line generated by your original 3D vector all you need to do is find its negative reciprocal, but that only works in 2D, and with that value you can use that and any arbitrary point on your original vector to generate a valid arbitrary perpendicular vector.

If you know your trigonometric functions and the properties of them; specific functions will have a value of 1, -1, 0, 1/2 as well as other specific values as output with corresponding specific values as input such as:

sin A { 0PI, PI, 2PI, 0deg, 180deg, 360deg } = 0
sin A { PI/2, 90deg } = 1 
sin A { 3PI/2, 270deg } = -1 
cos A { 0PI, 2PI, 0deg, 360deg } = 1 
cos A { PI/2, 3PI/2, 90deg, 270deg } = 0 
cos A { PI, 180deg } = -1 
tan A { 0PI, PI, 2PI, 0deg, 180deg, 360deg } = 0 
tan A { PI/4, 45deg } = 1
tan A { PI/2, 3PI/2, 90deg, 270deg } = undefined.

And with these functions and knowing the relationships of them and knowing that the product of two slopes = -1. We can see that the easiest of these functions to work with is the cos A where A = PI or 180 degrees and knowing the identity functions of trig such as the Reciprocal, Quotient, Pythagorean, Cofunction, (Even/Odd) Identities along with the various formulas such as Sum & Difference, Double Angle, Power Reducing, Sum to Product, and Product to Sum Formulas as well as the Arc functions and Inverse functions of the trigonometric functions you should be able to easily find a perpendicular line with the cosine version of the dot and cross products through the Law of Cosines. Here is more math about this: dot:cross and here wiki

So now that we know some of the math involved and their properties and we have a vector A(x,y,z) we can use the following:

// Values of Cosine At Specific Angles
cos (PI) = -1, cos (PI/2) = 0, cos (PI) = 1

// Orthogonal Vectors
A dot B = 0
// Parallel Vectors - Multiples Of Each Other
A(1,2)
B(2,4)
C(4,8)
etc.
// Angle Between To Vectors where neither vector is a 0 vector
cos (theta) = (A dot B) / (magnitude(A) * magnitude(B))

These give us the foundation to construct a basis vector called the unit vector or a normalized vector of original vector and from that normalized vector it is quite easy to calculate an orthogonal vector from it knowing the few properties listed above. To normalize a vector to get its unit vector the following formula or function can be used:

C++ Version

// -----------------------------------------------------------------------
// normalize()
// Make The Length Of This Vector Equal To One
inline void Vector3::normalize() {

    float magnitude;
    magnitude = sqrt( x*x + y*y + z*z );
    if ( magnitude <= Math::ZERO ) {  // Math::ZERO = (float)1e-7;
        x = 0.0f;
        x = 0.0f;
        x = 0.0f;

        return;
    }

    magnitude = 1 / magnitude;
    x *= magnitude;
    y *= magnitude;
    z *= magnitude;

} // Normalize

// -----------------------------------------------------------------------
// getCosAngle()
// Returns The cos(Angle) Value Between This Vector And Vector V. This
// Is Less Expensive Than Using GetAngle
inline float Vector3::getCosAngle( const Vector3 &v3, const bool normalized ) {
    // a . b = |a||b|cos(angle)
    // -> cos-1((a.b)/(|a||b|))

    // Make Sure We Do Not Divide By Zero
    float magnitudeA = length();
    if ( magnitudeA <= Math::ZERO ) {
        // This (A) Is An Invalid Vector
        return 0;
    }

    float value = 0;

    if ( normalized ) {
        // v3 Is Already Normalized
        value = dot(v3)/magnitudeA;
    } else {
        float magnitudeB = v3.length();
        if ( magnitudeB <= Math::ZERO) {
            // B Is An Invalid Vector
            return 0;
        }

        value = dot(v3)/(magnitudeA * magnitudeB);
    }

    // Correct Value Due To Rounding Problem
    Math::constrain( -1.0f, 1.0f, value );

    return value;

} // getCosAngle

// -----------------------------------------------------------------------
// length()
// Return The Length Of This Vector
inline float Vector3::length() const {

    return sqrtf( x*x + y*y + z*z );

} // length

// -----------------------------------------------------------------------
// Constrain()
// Prevent Value From Going Outside The Min, Max Range.
template<class T>
inline void Math::constrain( T min, T max, T &value ) {

    if ( value < min ) {
        value = min;
        return;
    }

    if ( value > max ) {
        value = max;
    }

} // constrain

And once you have a normalized vector which the overall length of the vector = 1. Calculating cross and dot products against this vector is quite cheap and to find one orthogonal to it is easy and to calculate it would be as follows

Vector3 original( 2,3,5 );
Vector3 arbitrary( 7,8,6 ); // Any vector that intersects original.

Vector3 findOrthogonalVector( const Vector3& original, const Vector3& arbitrary ) {
    // Test Both Vectors To See If Either Are the 0 Vector If So Just Return the Zero Vector And Be Done
    if ( original.isZero() || arbitrary.isZero() ) {
        return Vector3( 0.0f, 0.0f, 0.0f );
    }

    Vector3 orthogonal(); // default constructor set to (0,0,0);
    original.normalize(); // Total Length or Magnitude of 1 but Gives Direction
    arbitrary.normalize();

    // This will not give you values in degrees; it will give you the value of the cosine at a given angle but not the angle itself for example cos(90) = -1
    // This will give a value of -1 and not 90.
    float cosAngle = original.getCosAngle( arbitrary, true ); // True since arbitrary is normalized       

    // Calculate Actual Angle From Trig Functions
    float angle = someTrigFuncToGiveAngle( cosAngle ); 
    // Clamp Angle between [0,180]
    angle = clampAngle( angle, 0.0f, 180.0f ); 

    float desiredAngle = 0;
    // If a dot b = 0; where dot product = (ax*bx + ay*by + az*bz) then 
    if ( angle == Math::ZERO || angle == 180.0f ) {
        // original & arbitrary are parallel (oops)
        // Quick Fix or Quick Hack - Change X component by 1 and recursively call this function only in this case
        arbitrary.x += 1;
        return findOrthogonalVector( original, arbitrary );

    } else if ( angle == 90.0f ) { // (if in degrees otherwise in radians)
        // Already Orthogonal 
        return arbitrary;

    } else if ( angle < 90.0f ) {
        desiredAngle = 90.0f - angle;  

    } else if ( angle > 90.0f && angle < 180.0f ) {
        desiredAngle = angle - 90.0f;

    } 
    // With this new angle we should be able to find the orthogonal vector      
    // Rotate The Arbitrary Vector by desiredAngle
    return orthogonal = arbitrary.rotate( desiredAngle );
}

The Following Methods In The Algorithm Above Are Not Shown:

• someTrigFuncToGiveAngle( ... );
• clampAngle( ... );
• vector.rotate( ... );

These can be left for the OP to do as an assignment.

So the basic algorithm is as follows:

• Take Original Vector and Any Arbitrary Vector neither can be a zero vector
• If Either are 0 Vectors Then Just Return A Zero Vector And Done With Function
• Normalize both vectors
• Calculate the Cosine Angle Between Both Vectors
• Get Actual Angle From Cosine Of Angle
• Clamp Angle Between [0,180]
• Check if Angle = 0; If True - Parallel change x component by one and call again passing in the updated Arbitrary vector
• Check if Angle = 90; If True - Already Orthogonal return Arbitrary
• If Angle is < 90; Calculate the Desired Angle by subtracting it from 90 then Rotate Arbitrary Angle by it and return the Orthogonal Or Desired Vector
• If Angle is > 90; Calculate the Desired Angle by subtracting 90 from it then Rotate Arbitrary Angle by it and return the Orthogonal Or Desired Vector