VMath, part 1: Vectors

For my first actual content post, I thought I’d dig in to what I’ve written on the math side of things, and explain how it all works and why I chose the things I chose.  For this, and all future programming posts, the language of choice will be C++.

First thing to note, is that a lot of programmers who find themselves without a provided math library (such as is the case when working with OpenGL like me, or Vulkan) will usually end up using GLM.  There’s nothing particularly wrong with GLM, and if you aren’t especially math-inclined, then I honestly suggest you use it too.  The reason that I don’t use GLM is that, by its own admission, it’s designed for “the 90% of your code that does 10% of the work”.  That’s all fine and dandy until you start writing that other 10%, and now you need to optimize.  Without understanding how the math really works, and how your math library is doing its job, you can’t make those vital optimizations to the most important parts of the code.  Still, if you aren’t math-inclined, you can feel free to skip any VMath posts and simply use GLM as a drop-in replacement; the interface is very similar, by design.

Lets start off simple, a class for a 2-dimensional vector:

#ifndef HVH_VMATH_VEC2_H
#define HVH_VMATH_VEC2_H
#include <math.h>

namespace vmath {
class vec2
{}
} // namespace vmath
#endif // HVH_VMATH_VEC2_H

There’s really not much worth explaining here, yet.  Header guards (#ifndef, #define, #endif) are standard because they work on every compiler.  Even if you use Visual Studio, don’t use ‘#pragma once’ because it makes your code less portable to other compilers.

We use ‘namespace vmath’ because a lot of our classes use very simple and common names, and this is the best way of avoiding naming conflicts with other libraries.  Specifically, we use the same naming conventions as GLSL; that’ll come up later as well.

Next, we’ll add the meat of the class.  Its most important job is to store a pair of numbers which can be used for a point in 2d space, a direction and magnitude in 2d space, a color, or a texture co-ordinate, so the most important part is the data.

public:
    vec2() = default;
    vec2(const vec2&) = default;
    constexpr vec2(float x, float y) : x(x), y(y) {}

    union
    {
        struct { float x, y; };
        struct { float r, g; };
        struct { float s, t; };
        float data[2];
    };

First of all, you’ll notice how the first two constructors are explicitly declared ‘default’  The first allows us to construct a vec2 with no initialization arguments, but without defining any of its members; they’ll be whatever the memory happened to have there before, just like when you declare a ‘float’ on the stack without initializing it.  This is important because it allows ‘vec2’ to be a POD type.  You can read more about POD types HERE; long story short, it lets the compiler know that this is a simple type, and the compiler will allow us to do a few things because of that fact, including declare them as a ‘constexpr’, or constant expression.  The copy constructor , when declared default as we have, will simply value-copy each field into the new object, which is exactly what we want to happen.  The third constructor, the first normal one, lets us construct a ‘vec2’ using a list of values; we declare it ‘constexpr’ so we can use ‘vec2’ in ‘constexpr’ functions and values.

Now, whats this about this ‘union’ thing?  Think of it this way; our ‘vec2’ is three different structs, one which consists of ‘x’ and ‘y’, one which consists of ‘r’ and ‘g’, and one which consists of ‘s’ and ‘t’.  ‘union’ makes these structs overlap each other in memory, and have the exact same memory address.  When you access “my_vec2.x”, you’re accessing ‘x’, ‘r’, and ‘s’ at the same time.  This lets us use ‘vec2’ for different things and give the programmer some freedom about how they want to access their data in the most descriptive way possible, at no runtime or compile time cost whatsoever.

We use the same accessors that are available in GLSL: xyzw, rgba, and stpq.  ‘data’ allows us to access elements of the vector using an integer index, which we’ll be using later.

We good?  Let’s move on to operator overloads, a very handy feature of C++.

    // Equality test operators.
    inline constexpr bool operator == (vec2 rhs) const
        { return ((x == rhs.x) && (y == rhs.y)); }
    inline constexpr bool operator != (vec2 rhs) const
        { return ((x != rhs.x) || (y != rhs.y)); }

    // Arithmetic operators (vector with vector)
    inline constexpr vec2 operator + (vec2 rhs) const
        { return { x+rhs.x, y+rhs.y }; }
    inline constexpr vec2 operator - (vec2 rhs) const
        { return { x-rhs.x, y-rhs.y }; }
    inline constexpr vec2 operator * (vec2 rhs) const
        { return { x*rhs.x, y*rhs.y }; }
    inline constexpr vec2 operator / (vec2 rhs) const
        { return { x / rhs.x, y / rhs.y }; }

    // Arithmetic operators (vector with float)
    inline constexpr vec2 operator + (float rhs) const
        { return { x+rhs, y+rhs }; }
    inline constexpr vec2 operator - (float rhs) const
        { return { x-rhs, y-rhs }; }
    inline constexpr vec2 operator * (float rhs) const
        { return { x*rhs, y*rhs }; }
    inline constexpr vec2 operator / (float rhs) const
        { return { x/rhs, y/rhs }; }

    // Arithmetic operator (float with vector)
    friend inline constexpr vec2 operator + (const float lhs, const vec2 rhs)
        { return { lhs+rhs.x, lhs+rhs.y }; }
    friend inline constexpr vec2 operator - (const float lhs, const vec2 rhs)
        { return { lhs-rhs.x, lhs-rhs.y }; }
    friend inline constexpr vec2 operator * (const float lhs, const vec2 rhs)
        { return { lhs*rhs.x, lhs*rhs.y }; }
    friend inline constexpr vec2 operator / (const float lhs, const vec2 rhs)
        { return { lhs/rhs.x, lhs/rhs.y }; }

    // Arithmetic-assignment operators (vector with vector)
    inline vec2& operator += (vec2 rhs)
        { return *this = *this + rhs; }
    inline vec2& operator -= (vec2 rhs)
        { return *this = *this - rhs; }
    inline vec2& operator *= (vec2 rhs)
        { return *this = *this * rhs; }
    inline vec2& operator /= (vec2 rhs)
        { return *this = *this / rhs; }

    // Arithmetic-assignment operators (vector with float)
    inline vec2& operator += (float rhs)
        { return *this = *this + rhs; }
    inline vec2& operator -= (float rhs)
        { return *this = *this - rhs; }
    inline vec2& operator *= (float rhs)
        { return *this = *this * rhs; }
    inline vec2& operator /= (float rhs)
        { return *this = *this / rhs; }

    // Negation operator (unary minus)
    inline constexpr vec2 operator - () const
        { return { -x, -y }; }

We could go over each of these, but long story short, it simply allows you to use ‘vec2’ in a mathematical expression.  In linear algebra terms, addition, subtraction, multiplication, and division are all defined between two vectors and between vectors and scalars (single values), so we define both.  We need to use the ‘friend’ keyword for the operators where floats are on the left hand side of the expression, because operator overloads by default assume that ‘this’ (the vector) is the left hand side.

Now lets move on to some practical functionality, with member functions.

Lets add a simple test to see whether our vetor is zero:

inline constexpr bool is_zero() const
    { return ((x == 0.0f) && (y == 0.0f)); }

A function to give us the sum of each vector component:

inline constexpr float sum() const
    { return x + y; }

A function to give us the magnitude of the vector (this one isn’t constexpr because ‘sqrtf’ isn’t constexpr):

inline float magnitude() const
    { return sqrtf((x*x) + (y*y)); }

Sometimes, such as when comparing distances, the squared magnitude is enough, and much faster to compute:

inline constexpr float magnitude_squared() const
    { return (x*x) + (y*y); }

And, of course, normalization is important:

inline vec2 normalized() const
    { return *this / magnitude(); }

inline void normalize()
    { *this = normalized() }

We define a few functions as ‘static’ because they don’t really operate on a single ‘vec2’ instance.

‘minimum’ and ‘maximum’ give us the lowest and highest values, respectively, along each axis between two vectors:

static inline vec2 minimum(vec2 lhs, vec2 rhs)
    { return { fminf(lhs.x, rhs.x), fminf(lhs.y, rhs.y) }; }

static inline vec2 maximum(vec2 lhs, vec2 rhs)
    { return { fmaxf(lhs.x, rhs.x), fmaxf(lhs.y, lhs.y) }; }

We can use ‘clamp’ to clamp a vector’s values between two extremes:

inline vec2 clamped(vec2 lower, vec2 upper) const
    { return minimum(upper, maximum(lower, *this)); }

inline void clamp(vec2 lower, vec2 upper)
    { *this = clamped(lower, upper); }

Clamping the magnitude is also important; instead of clamping along each axis, we may want to ensure that our magnitude doesn’t exceed some distance.

inline vec2 clamped_magnitude(float max_magnitude) const
{
    float mag = magnitude();
    if (mag <= max_magnitude)
       return *this;
    else
       return this->normalized() * max_magnitude;
}

inline void clamp_magnitude(float max_magnitude)
    { *this = clamped_magnitude(max_magnitude); }

Knowing the distance between two points is helpful:

static inline float distance(vec2 a, vec2 b)
    { return (a-b).magnitude(); }

static inline float distance_squared(vec2 a, vec2 b)
    { return (a-b).magnitude_squared(); }

Moving a point towards another point with a maximum distance is useful:

inline vec2 moved_towards(vec2 target, float max_dist) const
{
    if (distance(*this, target) <= max_dist)
        return target;
    else
        return *this + ((target - *this).normalized() * max_dist);
}

inline void move_towards(vec2 target, float max_dist)
    { *this = moved_towards(target, max_dist); }

Find the dot product between two vectors:

static inline constexpr float dot(vec2 a, vec2 b)
        { return (a.x * b.x) + (a.y * b.y); }

Linearly interpolate between two vectors:

static inline vec2 lerp(vec2 from, vec2 to, float amount)
    { return (from * (1.0f - amount)) + (to * amount); }

And find the angle between two vectors:

static inline constexpr float determinant(vec2 a, vec2 b)
    { return (a.x * b.y) - (a.y * b.x); }

static inline float angle(vec2 from, vec2 to)
    { return atan2f(determinant(from, to), dot(from, to));

Frankly, a lot of these functions exist mostly for convenience.  They aren’t super optimized.  If I later go back and fix a bunch of them, I’ll make another post explaining how and why I changed it up.

Outside of the class, we can define a couple of constant vectors that we’ll use:

constexpr vec2 VEC2_ZERO = { 0, 0 };
constexpr vec2 VEC2_ONE = { 1, 1 };

Now, there’s one big feature that I developed for my vector class that I wanted to leave to the end, because it’s rather complicated, but it’s very handy to have.  That feature is swizzling.  If you’re familiar with GLSL, you’ve probably encountered it before; you can read about how it works HERE.  C++ was not designed with swizzling in mind, but we can fake it with the magic of templates.

The following code goes between the vector’s constructors and data definition:

private:
    friend class vec3;
    friend class vec4;

    // Here we define swizzle proxies.
    template <typename T, int X, int Y>
    struct SWIZZLE_READ_ONLY
    {
        inline operator vec2() const
        {
            T* self = (T*)this;
            return{ self->data[X], self->data[Y] };
        }

        inline vec2 operator + (vec2 rhs)
            { T* self = (T*)this; return{ self->data[X] + rhs.data[0], self->data[Y] + rhs.data[1] }; }
        inline vec2 operator - (vec2 rhs)
            { T* self = (T*)this; return{ self->data[X] - rhs.data[0], self->data[Y] - rhs.data[1] }; }
        inline vec2 operator * (vec2 rhs)
            { T* self = (T*)this; return{ self->data[X] * rhs.data[0], self->data[Y] * rhs.data[1] }; }
        inline vec2 operator / (vec2 rhs)
            { T* self = (T*)this; return{ self->data[X] / rhs.data[0], self->data[Y] / rhs.data[1] }; }

        inline vec2 operator + (float rhs)
            { T* self = (T*)this; return{ self->data[X] + rhs, self->data[Y] + rhs }; }
        inline vec2 operator - (float rhs)
            { T* self = (T*)this; return{ self->data[X] - rhs, self->data[Y] - rhs }; }
        inline vec2 operator * (float rhs)
            { T* self = (T*)this; return{ self->data[X] * rhs, self->data[Y] * rhs }; }
        inline vec2 operator / (float rhs)
            { T* self = (T*)this; return{ self->data[X] / rhs, self->data[Y] / rhs }; }
    };

    template <typename T, int X, int Y>
    struct SWIZZLE_READWRITE : public SWIZZLE_READ_ONLY<T, X, Y>
    {
        static_assert(X != Y, "Read+Write swizzling cannot be defined for multiple identical elements.");

        inline T& operator = (vec2 rhs)
        {
            T* self = (T*)this;
            self->data[X] = rhs.data[0];
            self->data[Y] = rhs.data[1];
            return *self;
        }
        
        // Any non-const, non-static member functions need to be defined here.

        inline vec2& operator += (vec2 rhs)
            { T* self = (T*)this; self->data[X] += rhs.data[0]; self->data[Y] += rhs.data[1]; return (vec2&)*self; }
        inline vec2& operator -= (vec2 rhs)
            { T* self = (T*)this; self->data[X] -= rhs.data[0]; self->data[Y] -= rhs.data[1]; return (vec2&)*self; }
        inline vec2& operator *= (vec2 rhs)
            { T* self = (T*)this; self->data[X] *= rhs.data[0]; self->data[Y] *= rhs.data[1]; return (vec2&)*self; }
        inline vec2& operator /= (vec2 rhs)
            { T* self = (T*)this; self->data[X] /= rhs.data[0]; self->data[Y] /= rhs.data[1]; return (vec2&)*self; }

        inline vec2& operator += (float rhs)
            { T* self = (T*)this; self->data[X] += rhs; self->data[Y] += rhs; return (vec2&)*self; }
        inline vec2& operator -= (float rhs)
            { T* self = (T*)this; self->data[X] -= rhs; self->data[Y] -= rhs; return (vec2&)*self; }
        inline vec2& operator *= (float rhs)
            { T* self = (T*)this; self->data[X] *= rhs; self->data[Y] *= rhs; return (vec2&)*self; }
        inline vec2& operator /= (float rhs)
            { T* self = (T*)this; self->data[X] /= rhs; self->data[Y] /= rhs; return (vec2&)*self; }

        inline void normalize()
        {
            T* self = (T*)this;
            vec2 temp = { self->data[X], self->data[Y] };
            temp.normalize();
            self->data[X] = temp.data[0]; self->data[Y] = temp.data[1];
        }

        inline void clamp(vec2 minimum, vec2 maximum)
        {
            T* self = (T*)this;
            vec2 temp = { self->data[X], self->data[Y] };
            temp.clamp(minimum, maximum);
            self->data[X] = temp.data[0]; self->data[Y] = temp.data[1];
        }

        inline void clamp_magnitude(float max_magnitude)
        {
            T* self = (T*)this;
            vec2 temp = { self->data[X], self->data[Y] };
            temp.clamp_magnitude(max_magnitude);
            self->data[X] = temp.data[0]; self->data[Y] = temp.data[1];
        }

        inline void move_towards(vec2 target, float max_dist)
        {
            T* self = (T*)this;
            vec2 temp = { self->data[X], self->data[Y] };
            temp.move_towards(target, max_dist);
            self->data[X] = temp.data[0]; self->data[Y] = temp.data[1];
        }
    };

public:

We declare all this stuff ‘private’ because we don’t want other code to be able to access any of it.  We declare ‘vec3’ and ‘vec2’ as “friends” because we DO want THEM to be able to access it.  Next, we change the union where our data is stored to:

//// Data storage
union
{
    struct { float x, y; };
    struct { float r, g; };
    struct { float s, t; };
    float data[2];

    //// Swizzles
    // xy
    SWIZZLE_READ_ONLY<vec2, 0, 0> xx;
    SWIZZLE_READWRITE<vec2, 0, 1> xy;
    SWIZZLE_READ_ONLY<vec2, 1, 1> yy;
    SWIZZLE_READWRITE<vec2, 1, 0> yx;
    // rg
    SWIZZLE_READ_ONLY<vec2, 0, 0> rr;
    SWIZZLE_READWRITE<vec2, 0, 1> rg;
    SWIZZLE_READ_ONLY<vec2, 1, 1> gg;
    SWIZZLE_READWRITE<vec2, 1, 0> gr;
    // st
    SWIZZLE_READ_ONLY<vec2, 0, 0> ss;
    SWIZZLE_READWRITE<vec2, 0, 1> st;
    SWIZZLE_READ_ONLY<vec2, 1, 1> tt;
    SWIZZLE_READWRITE<vec2, 1, 0> ts;
};

Alright, well that all looks like a mess.  What exactly is happening here?  Basically, SWIZZLE_READ_ONLY and SWIZZLE_READWRITE are proxy classes, which take advantage of the fact that they exist in the same memory as the real vector.  When created, they’re given integers that tell them which element goes in which place. We have two different classes because a swizzle that has two of the same element can’t be written to; GLSL disallows this as well.

‘operator vec2()’ and ‘operator = ()’  should do most of the heavy lifting, converting the proxy into a real vector to be used wherever.  We also need to implement any non-const operators and functions in SWIZZLE_READWRITE so that the data in the original vector will be modified. Theoretically this should be enough, but Visual Studio at least will complain when you try to perform arithmetic without the re-implemented operators, so we implement those in SWIZZLE_READ_ONLY.  SWIZZLE_READWRITE is a child of SWIZZLE_READ_ONLY simply so we don’t have to write those functions twice.

Whew, and that’s everything!  One complete ‘vec2’ class!  Best part is, that’s really all there is to it.  For ‘vec3’ and ‘vec4’, it’s really just a matter of handling the extra dimensions exactly as they’re handled here.  There are a couple of noteworthy exceptions to that rule, of course.  ‘vec3’ has quite a bit of additional functionality, such as the cross product, and handles angles differently.  ‘vec4’ doesn’t handle angles at all.  For the sake of completeness, let’s go over the noteworthy exceptions here.

First, most notably, the cross product, which is unique to ‘vec3’.  This creates a vector perpendicular to the two input vectors, and is very useful when working in 3 dimensions:

// Cross product: A x B = (a2b3 - a3b2, a3b1 - a1b3, a1b2 - a2b1)
static inline constexpr vec3 cross(vec3 a, vec3 b)
    { return { (a.y*b.z) - (a.z*b.y), (a.z*b.x) - (a.x*b.z), (a.x*b.y) - (a.y*b.x) }; }

Next, ortho-normalize, which uses the cross-product to make the two input vectors be perpendicular to each other:

static inline vec3 ortho_normalize(vec3 a, vec3 b)
{
    vec3 r = cross(a.normalized(), b.normalized());
    return cross(r, a.normalized());
}

Angle, which returns the shortest positive angle between the two vectors:

static float angle(vec3 a, vec3 b)
    { return acosf(clampf(dot(a, b) / (a.magnitude() * b.magnitude()), -1.0f, 1.0f)); }

And signed angle, for when it’d be kinda nice to tell whether you’re swinging left or right:

static float signed_angle(vec3 a, vec3 b)
{
    float my_angle = angle(a, b);
    if (cross(a, b).y < 0) my_angle = -my_angle;
    return my_angle;
}

Finally, I added this ‘hsv_to_rgb’ function, because vec3’s are often used for colors and, well, this just seemed like the best place for it:

// Transforms hue/saturation/value to red/green/blue.
static inline vec3 hsv_to_rgb(float hue, float sat, float val)
{
    vec3 result;
    result.r = fabsf(hue * 6.0f - 3.0f) - 1.0f;
    result.g = 2 - fabsf(hue * 6.0f - 2.0f);
    result.b = 2 - fabsf(hue * 6.0f - 4.0f);
    result.clamp({ 0,0,0 }, { 1,1,1 });
    return ((result - 1.0f) * sat + 1.0f) * val;
    }

And, since vec3 is probably the most used vector we’ll encounter, we define a few more constants:

constexpr vec3 VEC3_ZERO = { 0, 0, 0 };
constexpr vec3 VEC3_ONE = { 1, 1, 1 };
constexpr vec3 VEC3_UP = { 0, 1, 0 };
constexpr vec3 VEC3_DOWN = { 0, -1, 0 };
constexpr vec3 VEC3_FORWARD = { 0, 0, -1 };
constexpr vec3 VEC3_BACK = { 0, 0, 1};
constexpr vec3 VEC3_RIGHT = { 1, 0, 0 };
constexpr vec3 VEC3_LEFT = { -1, 0, 0 };

constexpr vec3 VEC3_BLACK = { 0, 0, 0 };
constexpr vec3 VEC3_WHITE = { 1, 1, 1 };
constexpr vec3 VEC3_RED = { 1, 0, 0 };
constexpr vec3 VEC3_GREEN = { 0, 1, 0 };
constexpr vec3 VEC3_BLUE = { 0, 0, 1 };
constexpr vec3 VEC3_YELLOW = { 1, 1, 0 };
constexpr vec3 VEC3_CYAN = { 0, 1, 1 };
constexpr vec3 VEC3_MAGENTA = { 1, 0, 1 };

Despite the extra dimension, ‘vec4’ is actually a lot simpler.  It’s basically ‘vec2’ but without the ‘angle’ stuff.

I’ll go ahead and share the source files for you here, so you can examine how it works and maybe even use it in your own code.  This code is provided as-is, with no guarantees or warranty.  Public domain, I don’t care if you use it in commercial projects, yada yada.

vmath_includes.h (keeps all the standard headers and a handful of common functions in a single place)

vec2.h

vec3.h

vec4.h

There’s a handful of functions in there that I didn’t go over, namely the stuff with springs and smooth damping.  That has to do more with calculus than linear algebra; it’s kinda physics-ey.  Honestly, I don’t fully understand it yet myself, I was just playing around with it.  Once I understand it properly and have a chance to flesh it out some, I’ll make another post going into detail.  For now, it’s safe to ignore.

Next part wont be quite so long, I promise!

 

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s