Tag: c-sharp

Operator Overloading in C++

One of the advantages of working with C++ is the ability of the programmer to overload not only functions, but also operators as well.  With operator overloading most operators can be extended to work not only with built-in types like floats and ints but also classes.  This gives the programmer the freedom to change how a built-in operator is used on objects of that particular class.  For example, we can overload the plus, minus and star operators to perform vector addition, subtraction and scalar multiplication on the components of our Vector3 class.  (I will be using this Vector3 class as an example of operator overloading for the remainder of this post.  For more on vector math see Vectors Part 1.)

The following code sample is an example of a three  dimensional vector class which uses operator overloading for many the the various operators that are available.  We will examine the syntax of each of these operators individually.

class Vector3 {
public:
    float x;
    float y;
    float z;

    //Constructor method
    Vector3(float x, float y, float z) : x(x), y(y), z(z) {}

    //unary operations
    Vector3 operator- () const { return Vector3(-x, -y, -z); }

    //binary operations
    Vector3 operator- (const Vector3& rhs) const {
        return Vector3(x - rhs.x, y - rhs.y, z - rhs.z);
    }
    Vector3 operator+ (const Vector3& rhs) const {
        return Vector3(x + rhs.x, y + rhs.y, z + rhs.z);
    }
    Vector3 operator* (const float& rhs) const {
        return Vector3(x * rhs, y * rhs, z * rhs);
    }
    friend Vector3 operator* (const float& lhs, const Vector3& rhs) {
        return rhs * lhs;
    }
    Vector3 operator/ (const float& rhs) const {
        return Vector3(x / rhs, y / rhs, z / rhs);
    }
    bool operator!= (const Vector3& rhs) const {
        return (*this - rhs).sqrMagnitude() >= 1.0e-6;
    }
    bool operator== (const Vector3& rhs) const {
        return (*this - rhs).sqrMagnitude() < 1.0e-6;
    }

    //assignment operation
    Vector3& operator= (const Vector3& rhs)  {
        //Check for self-assignment
        if (this == &rhs)
            return *this;
        x = rhs.x;
        y = rhs.y;
        z = rhs.z;
        return *this;
    }

    //compound assignment operations
    Vector3& operator+= (const Vector3& rhs) {
        x += rhs.x;
        y += rhs.y;
        z += rhs.z;
        return *this;
    }
    Vector3& operator-= (const Vector3& rhs) {
        x -= rhs.x;
        y -= rhs.y;
        z -= rhs.z;
        return *this;
    }
    Vector3& operator*= (const float& rhs) {
        x *= rhs;
        y *= rhs;
        z *= rhs;
        return *this;
    }
    Vector3& operator/= (const float& rhs) {
        x /= rhs;
        y /= rhs;
        z /= rhs;
        return *this;
    }

    //subscript operation
    float& operator[] (const int& i) {
        if (i < 0 || i > 2) throw std::out_of_range("Out of Vector3 range\n");
        return (i == 0) ? x : (i == 1) ? y : z;
    }
    const float& operator[] (const int& i) const {
        if (i < 0 || i > 2) throw std::out_of_range("Out of Vector3 range\n");
        return (i == 0) ? x : (i == 1) ? y : z;
    }

    //typecast operations
    operator Vector2() { return Vector2(x, y); }
    operator Vector4() { return Vector4(x, y, z, 0.0f); }

    ...

};

Unary Operators:

First let us take a look at the negate (-) operator.  The body of this operator simply returns a new Vector3 object with each of its floating point components negated.  But what we are most interested in is its declaration.  Like a typical function it has a return type which is Vector3 and a parameter list enclosed by parenthesis.  But instead of a function name we follow the return type with the keyword operator and then the symbol for the particular operator that we wish to overload, which in this case is the minus sign.

Vector3 operator- () const { return Vector3(-x, -y, -z); }

Overloaded operators can be defined either as class member functions or as non-member functions.  Here we define them as class member functions.  This means that we do not have to include one of the objects involved in the operation in the parameter list since we can reference that object through the this operator  (we could have just as easily written the x parameter in the body of the function as this->x).  The parameter list for our unary operator is empty because we are only operating on the current object.  If we were to write this as a non-member function its declaration would look something like this:

Vector3 operator- (const Vector3& a) { return Vector3(-a.x, -a.y, -a.z); }

Note that the const written after the parameter list in the class member function simply ensures that the function cannot make changes to the original object.  The same functionality is accomplished by the presence of the const keyword in the parameter list of the non-member function.  If you do wish the operator to make changes to the original object then you would simply remove the const in either of these places.

Binary/Comparison Operators:

Next let us observe the binary operations of addition and subtraction.

Vector3 operator- (const Vector3& rhs) const {
    return Vector3(x - rhs.x, y - rhs.y, z - rhs.z);
}
Vector3 operator+ (const Vector3& rhs) const {
    return Vector3(x + rhs.x, y + rhs.y, z + rhs.z);
}

As you can see their declarations are hardly any different from the unary operator except that they have an additional parameter which we have named rhs to indicate that it is the object that appears on the right-hand side of the operator.  Note that the parameter is passed as a constant reference.  We already discussed what purpose the const holds in this context, but why pass as reference?  Often when we pass an object by reference it is because we want to be able to make changes to the object within the body of the function, but in this case we pass it as constant to prevent that from happening.  So why do we do this?  The answer is simple.  Passing an object by value means that we have to take extra time to make a copy of the object that is passed.  Passing by reference allows us to skip that step by simply passing the address of the object, thus making our code more efficient.  Comparison operators generally work the same way with the exception that they simply return a boolean value.

Now let us take a look at our multiplication operator.

Vector3 operator* (const float& rhs) const {
    return Vector3(x * rhs, y * rhs, z * rhs);
}
friend Vector3 operator* (const float& lhs, const Vector3& rhs) {
    return rhs * lhs;
}

The first thing we notice is that we actually have two multiplication operator functions.  The first one looks just like our addition/subtraction operator declarations except that the parameter is of type float.  The second one on the other hand is actually a non-member function that we declare to be a friend of the Vector3 class to ensure that it has access to the Vector3 class’ private methods and variables.  Although this class does not have any private methods or variables we still need to declare this function as a friend because otherwise the compiler expects the function to have only one parameter.  But why do we need both of these functions in the first place?  Strictly speaking, we don’t need both unless we wish to ensure that scalar multiplication is a commutative property of our Vector3 class.  Our first function does a scalar multiplication operation for the case where we have a Vector3 object on the left and a scalar (float) on the right of our operator but would not work if we were to reverse the order.  That is why we need the second function to perform the operation when the order is reversed.  Note that we did not need to do the same thing for the addition operator despite it also being commutative because both the left and the right hand sides are of type Vector3.

Assignment Operators:

Next we will see how we can override the assignment operator.

Vector3& operator= (const Vector3& rhs)  {
        //Check for self-assignment
        if (this == &rhs)
            return *this;
        x = rhs.x;
        y = rhs.y;
        z = rhs.z;
        return *this;
    }

Most of the time we do not really need to override the assignment operator because the compiler generated constructor and assignment operator are usually sufficient.  But if our class contains pointers the default assignment operator can sometimes lead to problems because we end up with objects that have pointers to the same location.  To avoid these problems you will want to know how to override the assignment operator so that it makes a deep copy as well as how to make a copy constructor.  I will not be covering how to do either of these things in this post.  The example that I give here has no practical purpose since the default assignment operator does exactly the same thing.  However, there are two things worth mentioning if you ever do decide to overload the assignment operator.  The first is that, if you want your class to support chain assignment (see example below), the function should return the address of the current object.  The second is that you should always check for self-assignment before altering any data otherwise your class may end up releasing the resources that it is trying to copy from.

Compound assignment operators are generally simpler to understand.  They combine normal binary operations with an assignment operator.

Vector3& operator+= (const Vector3& rhs) {
    x += rhs.x;
    y += rhs.y;
    z += rhs.z;
    return *this;
}

Unlike normal binary operators they are generally meant to make changes to the original object, but like the assignment operator they usually return the address of the object.  This allows the programmer to write statements like the following:

a = b += c -= d *= e;  //example of chain assignment

Subscript/Function Call Operators:

Next we have the subscript operator.

float& operator[] (const int& i) {
    if (i < 0 || i > 2) throw std::out_of_range("Out of Vector3 range\n");
    return (i == 0) ? x : (i == 1) ? y : z;
}
const float& operator[] (const int& i) const {
    if (i < 0 || i > 2) throw std::out_of_range("Out of Vector3 range\n");
    return (i == 0) ? x : (i == 1) ? y : z;
}

Here I define it so that the subscript operator returns one of the vector components based on the integer value passed.  Notice however that once again we have two functions.  The first behaves just as we would expect.  The second one on the other hand is written as a constant function that returns a constant reference.  Here the const at the end of the second member function allows the function to be called even on an object that is defined as constant.  Without this function the programmer would not be able to access the contents of a constant object using the subscript notation.  The reason we return a constant reference is to ensure that the caller cannot modify the contents of our constant vector object.

If we were overloading the function call operator we would similarly wish to define two functions, one for constant objects and one for non-constant objects.  We do not override the function call operator in the Vector3 class, but something to keep in mind is that, unlike other operators, the function call operator can have as many parameters as the programmer likes.  An example where this might be useful is in a matrix class where the programmer wishes to look up a matrix element by its row and column.  You cannot do this with the subscript operator because it only takes one parameter, but it can be done in the following way:

class Matrix4x4 {
private:
    float data[16];
public:
    float& operator() (const int& row, const int& col) {
        return data[row + (col * 4)];
    }
    const float& operator() (const int& row, const int& col) const {
        return data[row + (col * 4)];
    }

    ...
};

Typecast Operators:

Operator overloading can even be used to define typecasting operations for our user-defined classes.  Suppose elsewhere we have defined a two-dimensional and a four-dimensional vector class and we want our compiler to know how to cast our Vector3 object into one of these types.  We can accomplish this in the following way:

operator Vector2() { return Vector2(x, y); }
operator Vector4() { return Vector4(x, y, z, 0.0f); }

Here we simply write the operator keyword followed by the type we want to cast our object into followed by open and close parentheses.  We do not need to define a return type since C++ assumes that we will return the correct type.  With the help of these functions we should be able write lines of code like the following without our compiler spitting out an error:

Vector2 v2 = Vector2(1.0f, 2.0f);
Vector3 v3 = Vector3(3.0f, 4.0f, 5.0f);
Vector4 v4 = Vector4(6.0f, 7.0f, 8.0f, 9.0f);
v2 = v3; //compiler implicitly typecasts to Vector2
v4 = v3; //compiler implicitly typecasts to Vector3

End of Lesson:

Well that is it for my lesson on operator overloading.  I hope the things you have learned here will be of use to you in designing your C++ applications.  If you are interested in reading some more of my postings, you may follow the links below:

 

Game Programming: Math Libraries

Most game engines define their own custom set of math libraries for use in their games.  Unity, for example, has defined its own Mathf library for  performing common operations on floating point numbers.  You may wonder what is wrong with just using System.Math as apposed to UnityEngine.Mathf.  Some people use Mathf instead of Math because it saves them the trouble of having to do double to float conversions.  However, in most cases Math is considered faster than Mathf.  Standard math libraries are certainly convenient and will most often be more efficient than any math libraries that you roll yourself, so why go to the trouble of implementing a custom library.  The problem with standard math libraries is that they may not always have the kinds of functions that are most important to game programming.  Other times it may be better to forgo accuracy in favor of efficiency by replacing expensive operations such as cos, sin and ex with less expensive and less accurate approximations of these functions.  In this article I will discuss some of the most commonly used functions found in Unity’s Mathf library and how you may implement them yourself.

Approximately:

Due to the nature of floating point numbers it is very common for operations performed on them to result in some form of rounding error.  This can lead to problems when trying to compare the equality of floating point numbers.  For example we may end up comparing a number like 2.0 to a number like 1.999999999.  For all effects and purposes these numbers are practically the same but when we check for equality we get a false statement.  This is why we sometimes need to take into account some margin of error when comparing floats.  This is exactly what the Mathf.Approximately function does.  It can be implemented in the following way:

bool Approximately(float a, float b) {
    return Abs(a - b) < Epsilon;
}

Here we take the absolute value of the difference between the two input variables and check if it is less than some arbitrarily small value called Epsilon.  It should be noted however that although some standard math already have a predefined value for epsilon, it may actually be better for the programmer to define their own value for epsilon since the standard math library’s definition of epsilon may be too small to make any difference from straight equality comparison.  An alternative approach could also be to define an Approximately function that takes the tolerance as a parameter and just replace Epsilon with this tolerance value.

Clamp:

Clamp is such a useful tool that it is any wonder why it is not a common tool found in every math library, especially since it is so simple to write.  Perhaps it is just so simple that some programmers feel it is not necessary to include it in a standard math library.  Simply put, clamp ensures that its input value lies within a certain range of values.  If the number lies between the given range it simply returns the number itself, otherwise it returns the nearest endpoint.  An example where this might be useful is in the case where the programmer wants to ensure that a player character never goes past the edges of the screen.  All the programmer has to do is Clamp the player’s position to within the screen boundaries.  Clamp can be implemented in the following way:

float Clamp(float value, float min, float max) {
    return (value < min)? min : (value > max)? max : value;
}

This function is so useful that Unity has even included a Clamp01 function that specifically clamps a floating point value to between 0 and 1.

Linear Interpolation (Lerp):

Linear interpolation is the process by which we find a point somewhere between two input values.  We can visualize it as a line drawn between two points.  Using the equation for a line we can write the formula for the desired value y as y = a + (b - a)x where a and b are the start and end values respectively and x represents what fraction of the way we are between our two values.  Note that we only get a value between a and b when x lies within the range [0, 1]; any other value for x results in some value outside the interval [a, b].  We also see that y = a when x = 0 and y = b when x = 1.  However, to avoid rounding error caused by floating point arithmetic, it is best to rewrite the equation as y = a(1 - x) + bx (this ensures that y = b when x = 1).  Depending on how you want your Lerp function to behave for certain values of x you may wish to clamp it to between 0 and 1 so as to never overshoot the endpoints.  However, in some rare cases it may be better to leave the value unclamped.  The following code sample shows how we may implement the former:

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

For more help learning the proper use of linear interpolation, see How to Lerp like a Pro by Robert Utter.

Move Towards:

While linear interpolation certainly has its uses, it is sometimes not entirely practical since it requires the programmer to know what point we started at and what fraction of the way between the two points we are at a given time.  But suppose instead we want to procedurally step some value towards some target using some kind of displacement value.  We want a way to find our new value while making sure that we do not overshoot our target.  To do this we simply compare our expected displacement to the difference between our current value and our target value.  If our expected displacement is greater than the actual distance to our target then we just return our target value.  Here is an example:

float MoveTowards(float current, float target, float maxDelta) {
    float delta = target - current;
    return (maxDelta >= delta) ? target : current + maxDelta;
}

Note that in the implementation above a negative value of maxDelta will cause our output to procedurally step away from  our target value.

Smooth Damp:

While MoveTowards generally serves its intended purpose it can often result in game mechanics that look somewhat jerky since it forces something to come to a sudden stop once it reaches its goal.  If we want our game to look somewhat cleaner we need to apply smoothing to our calculations.  One way that we can apply smoothing is to apply some kind of S-curve to give us that smooth ease-in/ease-out motion.  One commonly used formula is x = 6t^5 - 15t^4 + 10t^3 which gives us such an S-curve and can be used with little computing cost.

C2 S-curve

Another technique to achieve smoothing is through exponential decay.  This is a kind of ease-out method of smoothing and can be accomplished by the statement x = x + (desiredX - x) * 0.1f * timeDelta where timeDelta is the time step between updates.  The advantage of this type of smoothing is that you can smooth toward a changing target and you do not need to know how much time has elapsed since the start.  However, the problem with this type of smoothing is that the initial motion is sudden.

Exponential Decay Diagram

The approach we will be using, however, will rely on the model for a critically damped spring.  This model fixes the problem we have with the exponential decay model by maintaining the ease-in property found in the S-curve.

Crtically Damped Spring Diagram

A spring will always try to maintain a certain length by exerting a force that pushes it towards that desired length.  Without damping this system will oscillate back and forth past the desired length forever.  We say the system is critically damped when it returns to equilibrium in the least possible amount of time.  We can write the equation for the force acted on the spring as F = -k(x - x_d) - bv where k is the spring constant, x_d - x is the displacement from our desired position, b is the damping coefficient and v is the velocity of the end of the spring that works to slow the spring down.  Using Newton’s second law of motion, F = ma, and writing the velocity and acceleration as the first and second derivatives of the position, we can obtain the ordinary differential equation (ODE):

m\frac{d^2x}{dt^2} + b\frac{dx}{dt} + kx = kx_d    (1)

The next step requires us to know how to solve for an equation of x.  The above ODE will have a solution of the form x(t) = x_p(t) + x_c(t) where x_p(t) is a particular solution of the equation and x_c(t) is the complementary solution of the associated homogeneous DE.  We see that the particular solution is of the form x_p(t) = A where A is a constant since the right side of Equation 1 is constant.  By taking the first and second derivatives of x_p(t) and plugging them into Equation 1 we see that A = x_d.

That was our particular solution; now let us find our complementary solution.  The homogeneous DE gives us the characteristic equation mr^2 + br + k = 0 which has the roots:

r = \frac{-b \pm \sqrt{b^2 - 4mk}}{2m}

It can be shown that the system will be critically damped when our equation has the real repeated roots r = -\sqrt{\frac{k}{m}} which only occurs when b^2 = 4mk.  We can use this to simplify Equation 1:

\frac{d^2x}{dt^2} + 2\omega\frac{dx}{dt} + \omega^2(x_d - x) = 0 where \omega = \sqrt{\frac{k}{m}}    (2)

In the above formula \omega is the spring’s natural frequency and represents the stiffness or strength of the spring.  But now we also have our complementary solution which is x_c(t) = (c_0 + c_1t)e^{-\omega t}.  So our general solution is:

x(t) = x_d + (c_0 + c_1t)e^{-\omega t}

Solving for the constants c_0 and c_1 using the initial conditions x(0) = x_0 and x'(0) = v_0 gives us the actual solution:

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

This formula now tells us everything we need to calculate the new position and by taking the derivative we can find the new velocity of the object after after this time step:

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

Now the only parameters we need for our smoothing function are our current position and velocity, our target position our elapsed time and finally our smoothness factor \omega.  But how do we figure out what we should use for our smoothness factor?  It would be much more intuitive to control our function based on how long we want it to take to reach equilibrium.  If we define this smooth time as being the expected time to reach the target when at maximum velocity, then we solve for it quite simply by setting the acceleration to zero and the maximum velocity to \frac{x_d - x_0}{t_{sm}} where t_{sm} is the smooth time and plugging these values into Equation 2.  We then see that \omega = \frac{2}{t_{sm}}.

One final note before we write our smoothing function is that Equation 3 and 4 contain e^x.  We talked before about how operations such as these are expensive so the following code shows a polynomial that can be used to approximate e^x.  While I do not know exactly how this polynomial is produced, I would hazard a guess and say that it is probably constructed by approximating a higher order Taylor series polynomial down to a third order polynomial.  Finally to make our function exactly like the SmoothDamp function of Unity’s Math library we also need to clamp the speed according to some defined max speed.

float SmoothDamp(float current, float target, float& currentVelocity, float smoothTime,
        float maxSpeed = Infinity, float deltaTime = Time.deltaTime) {
    smoothTime = Max(0.0001f, smoothTime);
    float omega = 2.0f / smoothTime;
    float x = omega * deltaTime;
    float exp = 1.0f / (1.0f + x + 0.48 * x * x + 0.235 * x * x * x);
    float deltaX = target - current;
    float maxDelta = maxSpeed * smoothTime;

    // ensure we do not exceed our max speed
    deltaX = Clamp(deltaX, -maxDelta, maxDelta);
    float temp = (currentVelocity + omega * deltaX) * deltaTime;
    float result = current - deltaX + (deltaX + temp) * exp;
    currentVelocity = (currentVelocity - omega * temp) * exp;

    // ensure that we do not overshoot our target
    if (target - current > 0.0f == result > target) {
        result = target;
        currentVelocity = 0.0f;
    }
    return result;
}

Resources:

Kirmse, Andrew.  (2004).  Game Programming Gems 4.  Charles River Media.  ISBN 9781-58450-295-1.

Zill, Dennis G.  (2009).  A First Course in Differential Equations with Modeling Applications.  Brooks Cole.  ISBN 978-0-495-10824-5

External Links:

Critically Damped Spring Smoothing

How to Lerp like a Pro

Unity Scripting API: Mathf

Unity Forum: Formula behind SmoothDamp

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

A-Star Pathfinder (Unity)

During a team project for a college course on game engines, I worked on creating a game called Icebox Anchorage using the Unity3D game engine.  During this time I built a pathfinding system that I intended to be used in our game.  Sadly it never made it into the final version of the game because it turned out to be unnecessary.  However, despite the fact that it was never utilized, it is still a fully functional system.

Using a set of nodes placed strategically about a scene, I formed a point-based grid that would act as the map that my pathfinding system would traverse.  Then, using the A* search algorithm, enemy AIs would then have the ability to locate and chase the player even through winding corridors and mazes by locating a path to the node closest to the player.  See a video demonstration below:

A* Algorithm:

The A* star algorithm is perhaps the most commonly used algorithm for finding the shortest path between two points on a map or grid.  It is amalgamation  of two other pathfinding algorithms: Dijkstra’s algorithm and Greedy Best-First-Search.  Dijkstra’s algorithm works by procedurally visiting vertices in the graph in search of a shortest path.  It maintains a two lists of vertices called the open set and the closed set.  The closed set contains the list of vertices that have already been visited starting with only the starting vertex or node.  As it does so it stores a value g(n) representing the total distance traveled from the start node.  Meanwhile, the open set contains all the vertices adjacent to those that have already been visited.  Dijkstra’s algorithm works by repeatedly removing vertices from the open set starting with those that have the smallest g(n) and adding them to the closed set until a path to the goal has been found or it has been determined that no such path exists.  The problem with Dijkstra’s algorithm is that often most of the paths that it checks are useless making Dijkstra’s algorithm computationally expensive.

Greedy Best-First-Search works in a similar way in that it maintains an open and closed set of vertices.  The difference is that instead of visiting the vertices closest to the starting point it visits the vertices closest to the goal vertex.  This is calculated according to some estimate h(n) called a heuristic.  This is generally much more efficient than Dijkstra’s algorithm until it runs into an obstacle.  Unfortunately, because it has to backtrack to find a way around obstacles, Greedy Best-First-Search often does not return a shortest path.

The A* algorithm takes the best of both these algorithms to find the shortest path in a way that ignores unlikely paths.  For each vertex in the graph it keeps track of not only of the travel distance g(n) to the node n, but also that nodes heuristic value h(n).  When choosing which vertices to visit first it goes by the vertices with the lowest cost value represented by f(n) = g(n) + h(n).  This way it balances the two in a way that allows it to find a shortest path without being overly expensive.

Unity Implementation:

Two nodes are considered adjacent to one another if there is a straight line path between the two nodes.  Each of the path nodes keeps a list of all the nodes that are adjacent to them.  A separate GameObject called AI Map acts as the parent of all the path nodes in the level and contains the script that contains the A* algorithm.  Here is how it was written:

public List<GameObject> FindPath(GameObject start, GameObject end) {
    // the set of tentative nodes to be evaluated
    Dictionary<GameObject, Node> open = new Dictionary<GameObject, Node>();

    // the set of nodes already evaluated
    Dictionary<GameObject, Node> closed = new Dictionary<GameObject, Node>();

    // map of the path taken to each of the node
    Dictionary<GameObject, Node> previous = new Dictionary<GameObject, Node>();

    Node startNode = new Node();  // g_value is given a default of infinity
    startNode.g_value = 0;
    startNode.f_value = Distance(start, end);
    open[start] = startNode;
    previous[start] = startNode;

    while (open.Count > 0) {
        // find the key value pair from the open set with the lowest f-score
        var current = open.Aggregate((a, b) =>
            (a.Value.f_value < b.Value.f_value)? a : b);

        // have we reached the goal?
        if (current.Key == end) {
            // reconstruct path and return;
            List<GameObject> path = new List<GameObject>();
            path.Add(end);

            Node scan = current.Value;
            // scan the list backwards
            while (scan  != null && scan.came_from != null) {
                path.Insert(0, scan.came_from);
                scan = previous[scan.came_from];
            }

            return path;
        }

        open.Remove(current.Key);
        closed[current.Key] = current.Value;

        PathNode node = current.Key.GetComponent<PathNode>();
        foreach (GameObject neighbor in node.adjacentNodes) {
            // ignore if the node has already been processed
            if (closed.ContainsKey(neighbor)) {
                continue;
            }

            // calculate the distance to the neighbor via the current node
            float tentative_g = current.Value.g_value + Distance(neighbor, current.Key);

            // Is this in the open set
            Node currentNeighbor;
            if (!open.TryGetValue(neighbor, out currentNeighbor)) {
                // if not create a new node with the maximum g-value
                // and add to the open set and map
                currentNeighbor = new Node();
                open[neighbor] = currentNeighbor;
                previous[neighbor] = currentNeighbor;
            }

            // Check if the new g-score is lower than the current one
            if (tentative_g < currentNeighbor.g_value) {
                // if so update the open set with this newly discovered shortest path
                currentNeighbor.came_from = current.Key;
                currentNeighbor.g_value = tentative_g;
                currentNeighbor.f_value = tentative_g + Distance(neighbor, end);
            }
        }
    }

    // no path found so return an empty list
    return new List<GameObject>();
}

The enemies contain a Seeker script that calls on the AI Map to find a path to its selected target.  The target object as well as the points on the returned path require a SeekerTarget script.  While the path nodes are considered to be static targets, the player is considered a moving target and in order to find a path to it, it is necessary to determine which path node it is closest to.  Path nodes store direction vectors to each of their adjacent nodes as well as the midpoint between every pair of of nodes.  This data is used to calculate the Voronoi region that each path node represents.  Using the sign of the dot product of two vectors it is possible to determine which node the seeker target is closest to.