Tag: orthonormal

Vectors Part 2: Programming in C++

In my last post I covered the basics of vector arithmetic.  In it I wrote about vector addition, subtraction, scalar multiplication, magnitude, normalization, dot products and cross products.  This time will we be learning about some more advanced uses of vectors in game programming.  Some of these operations you may recognize as methods of the Vector3 class of the Unity Game Engine.  Here I will be explaining how you may implement some of these methods for yourself.

I will be presenting the sample code in C++.  I chose this language largely because it supports the concept of operator overloading.  It is possible to adapt the code to be written in languages like Java by creating free or static member functions for operations like addition, subtraction and scalar multiplication.  But if you wish to know more about operator overloading in C++, you may visit my blog post on Operator Overloading.

Projection:

When we have two vectors \vec{a} and \vec{b} we can visualize the first vector \vec{a} as the sum of two vectors \vec{a}_1 and \vec{a}_2 where \vec{a}_1 is the vector component of \vec{a} parallel to vector \vec{b} and \vec{a}_2 is vector component of \vec{a} orthogonal to vector \vec{b}.  Because \vec{a}_1 is in the same (or opposite) direction as \vec{b} it is called the vector projection of \vec{a} onto \vec{b} or proj_{\vec{b}}\vec{a} and its magnitude is called the scalar projection of \vec{a} onto \vec{b} or comp_{\vec{b}}\vec{a}.  Similarly \vec{a}_2 is called the vector rejection of \vec{a} from \vec{b}.

Vector Projection Diagram

Sometimes in video games it is necessary to calculate one of these component vectors.  Fortunately, since \vec{a} = \vec{a}_1 + \vec{a}_2 we can easily calculate \vec{a}_2 in terms of \vec{a}_1.  So we just need to find \vec{a}_1.  We can find the length of \vec{a}_1 using the definition of \cos{\theta} where \theta is the angle between \vec{a} and \vec{b}.  We know that \cos{\theta} is adjacent over hypotenuse.  So we have:

\cos{\theta} = \frac{\|\vec{a}_1\|}{\|\vec{a}\|} \implies comp_{\vec{b}}\vec{a} = \|\vec{a}_1\| = \|\vec{a}\|\cos{\theta}

If we recall the geometric definition of the dot product we can simplify this equation to:

comp_{\vec{b}}\vec{a} = \|\vec{a}\| \|\hat{b}\|\cos{\theta} = \vec{a} \bullet \hat{b} = \frac{\vec{a} \bullet \vec{b}}{\|\vec{b}\|}

Now that we have its magnitude we can find \vec{a}_1 by multiplying this by the unit vector in the direction of \vec{b}.  So we get:

proj_{\vec{b}}\vec{a} = \vec{a}_1 = \|\vec{a}_1\|\hat{b} = (\frac{\vec{a} \bullet \vec{b}}{\|\vec{b}\|})(\frac{\vec{b}}{\|\vec{b}\|}) = \frac{\vec{a} \bullet \vec{b}}{\|\vec{b}\|^2}\vec{b} = \frac{\vec{a} \bullet \vec{b}}{\vec{b} \bullet \vec{b}}\vec{b}

Reflection:

Now that we know how to calculate the vector projection, we can now calculate the reflection of a vector.  This is useful in computer graphics for creating mirror effects.  To calculate the reflection of a vector we need only the vector \vec{v} that we are reflecting and the normal \hat{n} of the surface that our vector is reflecting off.  For simplicity we will assume that the normal vector is normalized.

Vector Reflection Diagram

We can see by the diagram that proj_{\hat{n}}\vec{r} is the same as proj_{\hat{n}}-\vec{v}.  Similarly, the projection of \vec{r} onto the plane whose normal is \hat{n} is the same as the projection of \vec{v} off the plane.  This projection can be found by simply subtracting proj_{\hat{n}}\vec{v} from \vec{v}.  So by summing the components together we get:

\vec{r} = proj_{\hat{n}}-\vec{v} + (\vec{v} - proj_{\hat{n}}\vec{v}) = (-\vec{v} \bullet \hat{n})\hat{n} + \vec{v} - (\vec{v} \bullet \hat{n})\hat{n} = \vec{v} - 2(\vec{v} \bullet \hat{n})\hat{n}

Orthonormalization:

To understand the process of orthonormalization, we first need to understand the concept of vector spaces and bases.  A vector space V is quite simply a set, whose elements are vectors, for which we have defined two operations:  vector addition and scalar multiplication.  One example of V would be the n-dimensional Euclidean space \mathbb{R}^n.  For the following operation we will be looking specifically at \mathbb{R}^3.  A basis of a vector space is a subset of vectors \vec{v}_1, \dots, \vec{v}_n in V which span the vector space and are linearly independent, that is \vec{v}_i \neq k\vec{v}_j for any scalar constant k when i \neq j.  We say that a basis spans a vector space so long as any vector \vec{v} in the vector space can be uniquely written as:

\vec{v} = a_1\vec{v}_1 + a_2\vec{v}_2 + \dots + a_n\vec{v}_n

for some set of scalar values a_1, a_2, \dots , a_n.  An orthonormal basis is simply a basis whose vectors are normalized and whose inner product is \langle \vec{v}_i, \vec{v}_j \rangle = 0 when i \neq j.  The inner product is just a generalization of the dot product as it is applied to a set of vectors.  So what \langle \vec{v}_i, \vec{v}_j \rangle = 0 means is that all the vectors in the set are mutually orthogonal to one another.  For example the standard basis of the Cartesian coordinate system defines a point in 3D space by:

\vec{A} = x\hat{i} + y\hat{j} + z\hat{k} = x \begin{pmatrix}1 \\ 0 \\ 0\end{pmatrix} + y \begin{pmatrix}0 \\ 1 \\ 0\end{pmatrix} + z \begin{pmatrix}0 \\ 0 \\ 1\end{pmatrix} = \begin{pmatrix}x \\ y \\ z\end{pmatrix}

where the set of unit vectors (\hat{i}, \hat{j}, \hat{k}) are an orthonormal basis representing the x, y and z coordinates of Cartesian space.

Orthonormalization is the process by which we take a vector space basis and transform it into an orthonormal basis.  The following code shows how this can be achieved in two and three dimensions using vector rejection (projection in the orthogonal direction).

//OrthoNormalize in two dimensions
void OrthoNormalize(Vector3& normal, Vector3& tangent)  {
    normal.Normalize();
    tangent = (tangent - Vector3::Project(tangent, normal)).normalized();
}

//OrthoNormalize in three dimensions
void OrthoNormalize(Vector3& normal, Vector3& tangent, Vector3& binormal) {
    OrthoNormalize(normal, tangent);
    binormal = binormal - Vector3::Project(binormal, normal);
    binormal = (binormal - Vector3::Project(binormal, tangent)).normalized();
}

Linear Interpolation (Lerp):

Sometimes in games we will need to be able to find a point that is some fraction of the way between two  points.  The formula to achieve this is quite simply as we would use the exact same formula that we would use for two scalar values (See Game Programming: Math Libraries).  The formula is c = a(1 - t) + bt.  The only difference is that instead of scalar values we are using vectors for a and b.  Like with scalar values it is often a good idea to clamp our t value to somewhere between 0 and 1 so as to never overshoot the endpoints.  A code sample might look like this:

Vector3 Lerp(Vector3 a, Vector3 b, float t) {
    t = Clamp01(t);
    return a * (1 - t) + b * t;
}

Spherical Linear Interpolation (Slerp):

Standard linear interpolation works nicely when we are interested in interpolating between two points, but what if our vectors are used to represent directions instead.  In this case it might be more practical to use spherical interpolation.  Spherical interpolation allows us to smoothly interpolate between two orientations as if we are moving along the surface of a circle or sphere.  The general equation for the spherical interpolation of vectors is defined as:

\vec{v}' = \frac{\sin (1 - t)\theta}{\sin \theta}\vec{v}_1 + \frac{\sin t\theta}{\sin \theta}\vec{v}_2

However, when we write our function we need to perform a check to see whether the two input directions are parallel to one another.  This case is special because we do not know about which axis we should rotate our direction vector.  We can handle this situation in one of two ways: we can either do nothing and return one of the two input vectors or we can simply plug our parameters into our Lerp function and return the result.  How we implement it depends on how we want this function to behave.  Play around with it and see which implementation works best for you.  The following code sample uses Lerp to handle the special case:

Vector3 Slerp(Vector3 from, Vector3 to, float t) {
    float dot = Dot(from.normalized(), to.normalized());
    if (Mathf.Approximately(abs(dot), 1.0f, 1.0e-6)) {
        // use linear interpolation
        return Lerp(from, to, t);
    }

    float theta = acos(dot);
    float temp = sin(theta);
    return from * (sin(theta * (1.0 - t)) / temp) + to * (sin(theta * t) / temp);
}

For more information on the proper use of linear interpolation you can read the following blog post titled How to Lerp like a Pro written by Robert Utter.

Move Towards:

While linear interpolation certainly has its uses, it is sometimes better to procedurally step a vector towards some target point as opposed to trying to guess what fraction of the way between two points an object should be at a given time.  In this case we want a way to find our new position while making sure that we do not overshoot our target.  We can do this in much the same way we do it for scalar values (See Game Programming: Math Libraries).  To do this we simply compare our expected displacement to the actual distance between our current position and our target.  If our expected move distance is greater than the actual distance to our target then we just return our target vector.  Here is an example implementation:

Vector3 MoveTowards(Vector3 current, Vector3 target, float maxDistanceDelta) {
    Vector3 delta = target - current;
    float mag = delta.magnitude();
    return (maxDistanceDelta >= mag) ?
        target : current + maxDistanceDelta * (delta / mag);
}

One nice thing about this function is we can also use it to move a position vector away from a certain point by simply inputting a negative value for the maxDistanceDelta parameter.

Rotate Towards:

Similarly we may find that we want to procedurally step a direction vector towards a certain orientation and magnitude.  The magnitude is easy enough to step towards since we can simply use a MoveTowards function that works for scalar values instead of vectors.  Unfortunately, to step the rotation is a bit more complicated.  While it is not to hard to handle the cases where we have overshot our target rotation, the rest of the time rotating our vector is best handled using quaternions:

Vector3 RotateTowards(Vector3 current, Vector3 target,
        float maxRadiansDelta, float maxMagnitudeDelta) {
    float targetMag = target.magnitude();
    float currentMag = current.magnitude();
    Vector3 targetNorm = target / targetMag;
    Vector3 currentNorm = current / currentMag;
    float newMagnitude = Mathf.MoveTowards(currentMag, targetMag, maxMagnitudeDelta);

    float dot = Dot(currentNorm, targetNorm);
    if (Mathf.Approximately(abs(dot), 1.0f, 1.0e-6)) {
        // only change the magnitude
        return currentNorm * newMagnitude;
    }

    //check if we overshoot our rotation
    float angle = acos(dot) - maxRadiansDelta;
    if (angle <= 0.0f) {
        return targetNorm * newMagnitude;
    } else if (angle >= Mathf.PI) {
        // if maxRadiansDelta is negative we may be rotating away from target
        return -targetNorm * newMagnitude;
    }

    Quaternion q = Quaternion::AngleAxis(maxRadiansDelta, Cross(current, target));
    Vector3 p = q * current;
    return p.normalized() * newMagnitude;
}

In the code above we construct a Quaternion  object using the angle and axis of the rotation we wish to apply to our vector.  Quaternions have the special property that when multiplied by a vector the result is a rotated direction vector.  So we simply multiply our vector object by our Quaternion object and then set its magnitude before returning our final direction vector.  For more on Quaternion rotation visit my blog pages at:

[not yet available]

or see the the blog post titled Understanding Quaternions by Jeremiah van Oosten.

Smooth Damping:

MoveTowards works to move a vector towards a desired goal and then come to a sudden stop when the target is reached.  This often results in some jerky and unexpected behavior, especially when our velocity is relatively large.  Sometimes we may want our object to decelerate as we get closer to our target position.  This is the purpose of SmoothDamp.  We will be using the formula for a critically-damped spring and its derivative as the model for this algorithm.  The formula is:

x(t) = x_d + ((x_0 - x_d) + (v_0 + \omega(x_0 - x_d))t)e^{-\omega t}

and its derivative is:

v(t) = (v_0 - (v_0 + \omega(x_0 - x_d))\omega t)e^{-\omega t}

In the formulas above  x_0 and v_0 are the initial position and velocity, x_d is our target position, t is our elapsed time and \omega is the frequency of the spring which can be represented as 2 divided by the smooth time.  To understand how we get these formulas you can read my posting on Game Programming: Math Libraries.

Now we just need to adapt this formula for smoothing vectors.

Vector3 SmoothDamp(Vector3 current, Vector3 target, Vector3& currentVelocity,
        float smoothTime, float maxSpeed, float deltaTime) {
    // check if we are already at target;
    if (current == target) return target;

    smoothTime = Mathf.Max(0.0001f, smoothTime);
    float omega = 2.0f / smoothTime;
    float x = omega * deltaTime;
    float exp = 1.0f / (1.0f + x + 0.48f * x * x + 0.235f * x * x * x);
    Vector3 delta = current - target;
    float mag = delta.magnitude;
    float maxDelta = maxSpeed * smoothTime;

    // ensure we do not exceed our max speed
    float deltaX = Mathf.Min(mag, maxDelta);
    delta = (delta * deltaX) / mag;

    Vector3 temp = (currentVelocity + omega * delta) * deltaTime;
    currentVelocity = (currentVelocity - omega * temp) * exp;
    Vector3 result = current - delta + (delta + temp) * exp;

    // ensure that we do not overshoot our target
    if ((target - current).sqrMagnitude <= (result - current).sqrMagnitude) {
        result = target;
        currentVelocity = Vector3::zero;
    }
    return result;
}

Summary:

We have now covered a preponderance of the things that we can do with vectors, but they are not the only concept from linear algebra important to game programming.  Next time we will be covering the basics of matrices and how they are applied to three dimensional transformations of game objects and setting up the projection plane of our game camera.

Resources:

Lay, David C.  (2005).  Linear Algebra and its Applications Third Updated Edition.  Addison Wesley.  ISBN 978-0-321-28713-4.

External Links:

Orthonormal Basis

Scalar and Vector Projections pdf