Appendix A. C++ Vector Library

A decent vector library is useful in any 3D graphics application but doubly so when you’re using OpenGL ES 2.0, since it doesn’t provide functions for matrix math on its own. Review “Vector Beautification with C++” for an overview of the library listed in this appendix. In brief, it’s template-based C++, and it’s composed solely of three simple header files. The library defines a set of types that are named after GLSL’s vector types.

Disclaimer Regarding Performance

This code herein is designed for readability and ease of use. It does not attempt to leverage the ARM’s special instructions for squeezing out maximum performance from your CPU.

First- and second-generation iPhones and iPod touches have dedicated hardware called the VFP unit, which can help with vector math. The iPhone 3GS has an extended instruction set for these operations called NEON. If you’re interested in using these special instruction sets, I suggest taking a look at the math portion of the excellent oolong library, available here:

The gcc compiler may be able to generate some vector-oriented ARM instructions when maximum optimizations are enabled, but there are no guarantees. In Xcode, select ProjectEdit Project Settings, go to the Build tab, select Show All Settings, and look for the Optimization Level option under the Code Generation heading. Don’t assume anything until you look at the generated assembly code.

Vector.hpp

#pragma once
#include <cmath>

const float Pi = 4 * std::atan(1.0f);
const float TwoPi = 2 * Pi;

template <typename T>
struct Vector2 {
Vector2() {}
Vector2(T x, T y) : x(x), y(y) {}
T Dot(const Vector2& v) const
{
return x * v.x + y * v.y;
}
Vector2 operator+(const Vector2& v) const
{
return Vector2(x + v.x, y + v.y);
}
Vector2 operator-(const Vector2& v) const
{
return Vector2(x - v.x, y - v.y);
}
void operator+=(const Vector2& v)
{
*this = Vector2(x + v.x, y + v.y);
}
void operator-=(const Vector2& v)
{
*this = Vector2(x - v.x, y - v.y);
}
Vector2 operator/(float s) const
{
return Vector2(x / s, y / s);
}
Vector2 operator*(float s) const
{
return Vector2(x * s, y * s);
}
void operator/=(float s)
{
*this = Vector2(x / s, y / s);
}
void operator*=(float s)
{
*this = Vector2(x * s, y * s);
}
void Normalize()
{
float s = 1.0f / Length();
x *= s;
y *= s;
}
Vector2 Normalized() const
{
Vector2 v = *this;
v.Normalize();
return v;
}
T LengthSquared() const
{
return x * x + y * y;
}
T Length() const
{
return sqrt(LengthSquared());
}
const T* Pointer() const
{
return &x;
}
operator Vector2<float>() const
{
return Vector2<float>(x, y);
}
bool operator==(const Vector2& v) const
{
return x == v.x && y == v.y;
}
Vector2 Lerp(float t, const Vector2& v) const
{
return Vector2(x * (1 - t) + v.x * t,
y * (1 - t) + v.y * t);
}
template <typename P>
P* Write(P* pData)
{
Vector2* pVector = (Vector2*) pData;
*pVector++ = *this;
return (P*) pVector;
}
T x;
T y;
};

template <typename T>
struct Vector3 {
Vector3() {}
Vector3(T x, T y, T z) : x(x), y(y), z(z) {}
T Length()
{
return std::sqrt(x * x + y * y + z * z);
}
void Normalize()
{
float s = 1.0f / Length();
x *= s;
y *= s;
z *= s;
}
Vector3 Normalized() const
{
Vector3 v = *this;
v.Normalize();
return v;
}
Vector3 Cross(const Vector3& v) const
{
return Vector3(y * v.z - z * v.y,
z * v.x - x * v.z,
x * v.y - y * v.x);
}
T Dot(const Vector3& v) const
{
return x * v.x + y * v.y + z * v.z;
}
Vector3 operator+(const Vector3& v) const
{
return Vector3(x + v.x, y + v.y,  z + v.z);
}
void operator+=(const Vector3& v)
{
x += v.x;
y += v.y;
z += v.z;
}
void operator-=(const Vector3& v)
{
x -= v.x;
y -= v.y;
z -= v.z;
}
void operator/=(T s)
{
x /= s;
y /= s;
z /= s;
}
Vector3 operator-(const Vector3& v) const
{
return Vector3(x - v.x, y - v.y,  z - v.z);
}
Vector3 operator-() const
{
return Vector3(-x, -y, -z);
}
Vector3 operator*(T s) const
{
return Vector3(x * s, y * s, z * s);
}
Vector3 operator/(T s) const
{
return Vector3(x / s, y / s, z / s);
}
bool operator==(const Vector3& v) const
{
return x == v.x && y == v.y && z == v.z;
}
Vector3 Lerp(float t, const Vector3& v) const
{
return Vector3(x * (1 - t) + v.x * t,
y * (1 - t) + v.y * t,
z * (1 - t) + v.z * t);
}
const T* Pointer() const
{
return &x;
}
template <typename P>
P* Write(P* pData)
{
Vector3<T>* pVector = (Vector3<T>*) pData;
*pVector++ = *this;
return (P*) pVector;
}
T x;
T y;
T z;
};

template <typename T>
struct Vector4 {
Vector4() {}
Vector4(T x, T y, T z, T w) : x(x), y(y), z(z), w(w) {}
Vector4(const Vector3<T>& v, T w) : x(v.x), y(v.y), z(v.z), w(w) {}
T Dot(const Vector4& v) const
{
return x * v.x + y * v.y + z * v.z + w * v.w;
}
Vector4 Lerp(float t, const Vector4& v) const
{
return Vector4(x * (1 - t) + v.x * t,
y * (1 - t) + v.y * t,
z * (1 - t) + v.z * t,
w * (1 - t) + v.w * t);
}
const T* Pointer() const
{
return &x;
}
T x;
T y;
T z;
T w;
};

typedef Vector2<bool> bvec2;

typedef Vector2<int> ivec2;
typedef Vector3<int> ivec3;
typedef Vector4<int> ivec4;

typedef Vector2<float> vec2;
typedef Vector3<float> vec3;
typedef Vector4<float> vec4;


Matrix.hpp

#pragma once
#include "Vector.hpp"

template <typename T>
struct Matrix2 {
Matrix2()
{
x.x = 1; x.y = 0;
y.x = 0; y.y = 1;
}
Matrix2(const T* m)
{
x.x = m[0]; x.y = m[1];
y.x = m[2]; y.y = m[3];
}
vec2 x;
vec2 y;
};

template <typename T>
struct Matrix3 {
Matrix3()
{
x.x = 1; x.y = 0; x.z = 0;
y.x = 0; y.y = 1; y.z = 0;
z.x = 0; z.y = 0; z.z = 1;
}
Matrix3(const T* m)
{
x.x = m[0]; x.y = m[1]; x.z = m[2];
y.x = m[3]; y.y = m[4]; y.z = m[5];
z.x = m[6]; z.y = m[7]; z.z = m[8];
}
Matrix3(vec3 x, vec3 y, vec3 z) : x(x), y(y), z(z)
{
}
Matrix3 Transposed() const
{
Matrix3 m;
m.x.x = x.x; m.x.y = y.x; m.x.z = z.x;
m.y.x = x.y; m.y.y = y.y; m.y.z = z.y;
m.z.x = x.z; m.z.y = y.z; m.z.z = z.z;
return m;
}
const T* Pointer() const
{
return &x.x;
}
vec3 x;
vec3 y;
vec3 z;
};

template <typename T>
struct Matrix4 {
Matrix4()
{
x.x = 1; x.y = 0; x.z = 0; x.w = 0;
y.x = 0; y.y = 1; y.z = 0; y.w = 0;
z.x = 0; z.y = 0; z.z = 1; z.w = 0;
w.x = 0; w.y = 0; w.z = 0; w.w = 1;
}
Matrix4(const Matrix3<T>& m)
{
x.x = m.x.x; x.y = m.x.y; x.z = m.x.z; x.w = 0;
y.x = m.y.x; y.y = m.y.y; y.z = m.y.z; y.w = 0;
z.x = m.z.x; z.y = m.z.y; z.z = m.z.z; z.w = 0;
w.x = 0; w.y = 0; w.z = 0; w.w = 1;
}
Matrix4(const T* m)
{
x.x = m[0];  x.y = m[1];  x.z = m[2];  x.w = m[3];
y.x = m[4];  y.y = m[5];  y.z = m[6];  y.w = m[7];
z.x = m[8];  z.y = m[9];  z.z = m[10]; z.w = m[11];
w.x = m[12]; w.y = m[13]; w.z = m[14]; w.w = m[15];
}
Matrix4 operator * (const Matrix4& b) const
{
Matrix4 m;
m.x.x = x.x * b.x.x + x.y * b.y.x + x.z * b.z.x + x.w * b.w.x;
m.x.y = x.x * b.x.y + x.y * b.y.y + x.z * b.z.y + x.w * b.w.y;
m.x.z = x.x * b.x.z + x.y * b.y.z + x.z * b.z.z + x.w * b.w.z;
m.x.w = x.x * b.x.w + x.y * b.y.w + x.z * b.z.w + x.w * b.w.w;
m.y.x = y.x * b.x.x + y.y * b.y.x + y.z * b.z.x + y.w * b.w.x;
m.y.y = y.x * b.x.y + y.y * b.y.y + y.z * b.z.y + y.w * b.w.y;
m.y.z = y.x * b.x.z + y.y * b.y.z + y.z * b.z.z + y.w * b.w.z;
m.y.w = y.x * b.x.w + y.y * b.y.w + y.z * b.z.w + y.w * b.w.w;
m.z.x = z.x * b.x.x + z.y * b.y.x + z.z * b.z.x + z.w * b.w.x;
m.z.y = z.x * b.x.y + z.y * b.y.y + z.z * b.z.y + z.w * b.w.y;
m.z.z = z.x * b.x.z + z.y * b.y.z + z.z * b.z.z + z.w * b.w.z;
m.z.w = z.x * b.x.w + z.y * b.y.w + z.z * b.z.w + z.w * b.w.w;
m.w.x = w.x * b.x.x + w.y * b.y.x + w.z * b.z.x + w.w * b.w.x;
m.w.y = w.x * b.x.y + w.y * b.y.y + w.z * b.z.y + w.w * b.w.y;
m.w.z = w.x * b.x.z + w.y * b.y.z + w.z * b.z.z + w.w * b.w.z;
m.w.w = w.x * b.x.w + w.y * b.y.w + w.z * b.z.w + w.w * b.w.w;
return m;
}
Vector4<T> operator * (const Vector4<T>& b) const
{
Vector4<T> v;
v.x = x.x * b.x + x.y * b.y + x.z * b.z + x.w * b.w;
v.y = y.x * b.x + y.y * b.y + y.z * b.z + y.w * b.w;
v.z = z.x * b.x + z.y * b.y + z.z * b.z + z.w * b.w;
v.w = w.x * b.x + w.y * b.y + w.z * b.z + w.w * b.w;
return v;
}
Matrix4& operator *= (const Matrix4& b)
{
Matrix4 m = *this * b;
return (*this = m);
}
Matrix4 Transposed() const
{
Matrix4 m;
m.x.x = x.x; m.x.y = y.x; m.x.z = z.x; m.x.w = w.x;
m.y.x = x.y; m.y.y = y.y; m.y.z = z.y; m.y.w = w.y;
m.z.x = x.z; m.z.y = y.z; m.z.z = z.z; m.z.w = w.z;
m.w.x = x.w; m.w.y = y.w; m.w.z = z.w; m.w.w = w.w;
return m;
}
Matrix3<T> ToMat3() const
{
Matrix3<T> m;
m.x.x = x.x; m.y.x = y.x; m.z.x = z.x;
m.x.y = x.y; m.y.y = y.y; m.z.y = z.y;
m.x.z = x.z; m.y.z = y.z; m.z.z = z.z;
return m;
}
const T* Pointer() const
{
return &x.x;
}
static Matrix4<T> Identity()
{
return Matrix4();
}
static Matrix4<T> Translate(const Vector3<T>& v)
{
Matrix4 m;
m.x.x = 1; m.x.y = 0; m.x.z = 0; m.x.w = 0;
m.y.x = 0; m.y.y = 1; m.y.z = 0; m.y.w = 0;
m.z.x = 0; m.z.y = 0; m.z.z = 1; m.z.w = 0;
m.w.x = v.x; m.w.y = v.y; m.w.z = v.z; m.w.w = 1;
return m;
}
static Matrix4<T> Translate(T x, T y, T z)
{
Matrix4 m;
m.x.x = 1; m.x.y = 0; m.x.z = 0; m.x.w = 0;
m.y.x = 0; m.y.y = 1; m.y.z = 0; m.y.w = 0;
m.z.x = 0; m.z.y = 0; m.z.z = 1; m.z.w = 0;
m.w.x = x; m.w.y = y; m.w.z = z; m.w.w = 1;
return m;
}
static Matrix4<T> Scale(T s)
{
Matrix4 m;
m.x.x = s; m.x.y = 0; m.x.z = 0; m.x.w = 0;
m.y.x = 0; m.y.y = s; m.y.z = 0; m.y.w = 0;
m.z.x = 0; m.z.y = 0; m.z.z = s; m.z.w = 0;
m.w.x = 0; m.w.y = 0; m.w.z = 0; m.w.w = 1;
return m;
}
static Matrix4<T> Scale(T x, T y, T z)
{
Matrix4 m;
m.x.x = x; m.x.y = 0; m.x.z = 0; m.x.w = 0;
m.y.x = 0; m.y.y = y; m.y.z = 0; m.y.w = 0;
m.z.x = 0; m.z.y = 0; m.z.z = z; m.z.w = 0;
m.w.x = 0; m.w.y = 0; m.w.z = 0; m.w.w = 1;
return m;
}
static Matrix4<T> Rotate(T degrees)
{
T radians = degrees * 3.14159f / 180.0f;

Matrix4 m = Identity();
m.x.x =  c; m.x.y = s;
m.y.x = -s; m.y.y = c;
return m;
}
static Matrix4<T> Rotate(T degrees, const vec3& axis)
{
T radians = degrees * 3.14159f / 180.0f;

Matrix4 m = Identity();
m.x.x = c + (1 - c) * axis.x * axis.x;
m.x.y = (1 - c) * axis.x * axis.y - axis.z * s;
m.x.z = (1 - c) * axis.x * axis.z + axis.y * s;
m.y.x = (1 - c) * axis.x * axis.y + axis.z * s;
m.y.y = c + (1 - c) * axis.y * axis.y;
m.y.z = (1 - c) * axis.y * axis.z - axis.x * s;
m.z.x = (1 - c) * axis.x * axis.z - axis.y * s;
m.z.y = (1 - c) * axis.y * axis.z + axis.x * s;
m.z.z = c + (1 - c) * axis.z * axis.z;
return m;
}
static Matrix4<T> Ortho(T left, T right, T bottom, T top, T near, T far)
{
T a = 2.0f / (right - left);
T b = 2.0f / (top - bottom);
T c = -2.0f / (far - near);
T tx = (right + left) / (right - left);
T ty = (top + bottom) / (top - bottom);
T tz = (far + near) / (far - near);
Matrix4 m;
m.x.x = a; m.x.y = 0; m.x.z = 0; m.x.w = tx;
m.y.x = 0; m.y.y = b; m.y.z = 0; m.y.w = ty;
m.z.x = 0; m.z.y = 0; m.z.z = c; m.z.w = tz;
m.w.x = 0; m.w.y = 0; m.w.z = 0; m.w.w = 1;
return m;
}
static Matrix4<T> Frustum(T left, T right, T bottom, T top, T near, T far)
{
T a = 2 * near / (right - left);
T b = 2 * near / (top - bottom);
T c = (right + left) / (right - left);
T d = (top + bottom) / (top - bottom);
T e = - (far + near) / (far - near);
T f = -2 * far * near / (far - near);
Matrix4 m;
m.x.x = a; m.x.y = 0; m.x.z = 0; m.x.w = 0;
m.y.x = 0; m.y.y = b; m.y.z = 0; m.y.w = 0;
m.z.x = c; m.z.y = d; m.z.z = e; m.z.w = -1;
m.w.x = 0; m.w.y = 0; m.w.z = f; m.w.w = 1;
return m;
}
static Matrix4<T> LookAt(const Vector3<T>& eye,
const Vector3<T>& target,
const Vector3<T>& up)
{
Vector3<T> z = (eye - target).Normalized();
Vector3<T> x = up.Cross(z).Normalized();
Vector3<T> y = z.Cross(x).Normalized();

Matrix4<T> m;
m.x = Vector4<T>(x, 0);
m.y = Vector4<T>(y, 0);
m.z = Vector4<T>(z, 0);
m.w = Vector4<T>(0, 0, 0, 1);

Vector4<T> eyePrime = m * Vector4<T>(-eye, 1);
m = m.Transposed();
m.w = eyePrime;

return m;
}
vec4 x;
vec4 y;
vec4 z;
vec4 w;
};

typedef Matrix2<float> mat2;
typedef Matrix3<float> mat3;
typedef Matrix4<float> mat4;


Quaternion.hpp

#pragma once
#include "Matrix.hpp"

template <typename T>
struct QuaternionT {
T x;
T y;
T z;
T w;

QuaternionT();
QuaternionT(T x, T y, T z, T w);

QuaternionT<T> Slerp(T mu, const QuaternionT<T>& q) const;
QuaternionT<T> Rotated(const QuaternionT<T>& b) const;
QuaternionT<T> Scaled(T scale) const;
T Dot(const QuaternionT<T>& q) const;
Matrix3<T> ToMatrix() const;
Vector4<T> ToVector() const;
QuaternionT<T> operator-(const QuaternionT<T>& q) const;
QuaternionT<T> operator+(const QuaternionT<T>& q) const;
bool operator==(const QuaternionT<T>& q) const;
bool operator!=(const QuaternionT<T>& q) const;

void Normalize();
void Rotate(const QuaternionT<T>& q);

static QuaternionT<T> CreateFromVectors(const Vector3<T>& v0, const Vector3<T>& v1);
static QuaternionT<T> CreateFromAxisAngle(const Vector3<T>& axis, float radians);
};

template <typename T>
inline QuaternionT<T>::QuaternionT() : x(0), y(0), z(0), w(1)
{
}

template <typename T>
inline QuaternionT<T>::QuaternionT(T x, T y, T z, T w) : x(x), y(y), z(z), w(w)
{
}

// Ken Shoemake's famous method.
template <typename T>
inline QuaternionT<T> QuaternionT<T>::Slerp(T t, const QuaternionT<T>& v1) const
{
const T epsilon = 0.0005f;
T dot = Dot(v1);

if (dot > 1 - epsilon) {
QuaternionT<T> result = v1 + (*this - v1).Scaled(t);
result.Normalize();
return result;
}

if (dot < 0)
dot = 0;

if (dot > 1)
dot = 1;

T theta0 = std::acos(dot);
T theta = theta0 * t;

QuaternionT<T> v2 = (v1 - Scaled(dot));
v2.Normalize();

QuaternionT<T> q = Scaled(std::cos(theta)) + v2.Scaled(std::sin(theta));
q.Normalize();
return q;
}

template <typename T>
inline QuaternionT<T> QuaternionT<T>::Rotated(const QuaternionT<T>& b) const
{
QuaternionT<T> q;
q.w = w * b.w - x * b.x - y * b.y - z * b.z;
q.x = w * b.x + x * b.w + y * b.z - z * b.y;
q.y = w * b.y + y * b.w + z * b.x - x * b.z;
q.z = w * b.z + z * b.w + x * b.y - y * b.x;
q.Normalize();
return q;
}

template <typename T>
inline QuaternionT<T> QuaternionT<T>::Scaled(T s) const
{
return QuaternionT<T>(x * s, y * s, z * s, w * s);
}

template <typename T>
inline T QuaternionT<T>::Dot(const QuaternionT<T>& q) const
{
return x * q.x + y * q.y + z * q.z + w * q.w;
}

template <typename T>
inline Matrix3<T> QuaternionT<T>::ToMatrix() const
{
const T s = 2;
T xs, ys, zs;
T wx, wy, wz;
T xx, xy, xz;
T yy, yz, zz;
xs = x * s;  ys = y * s;  zs = z * s;
wx = w * xs; wy = w * ys; wz = w * zs;
xx = x * xs; xy = x * ys; xz = x * zs;
yy = y * ys; yz = y * zs; zz = z * zs;
Matrix3<T> m;
m.x.x = 1 - (yy + zz); m.y.x = xy - wz;  m.z.x = xz + wy;
m.x.y = xy + wz; m.y.y = 1 - (xx + zz); m.z.y = yz - wx;
m.x.z = xz - wy; m.y.z = yz + wx;  m.z.z = 1 - (xx + yy);
return m;
}

template <typename T>
inline Vector4<T> QuaternionT<T>::ToVector() const
{
return Vector4<T>(x, y, z, w);
}

template <typename T>
QuaternionT<T> QuaternionT<T>::operator-(const QuaternionT<T>& q) const
{
return QuaternionT<T>(x - q.x, y - q.y, z - q.z, w - q.w);
}

template <typename T>
QuaternionT<T> QuaternionT<T>::operator+(const QuaternionT<T>& q) const
{
return QuaternionT<T>(x + q.x, y + q.y, z + q.z, w + q.w);
}

template <typename T>
bool QuaternionT<T>::operator==(const QuaternionT<T>& q) const
{
return x == q.x && y == q.y && z == q.z && w == q.w;
}

template <typename T>
bool QuaternionT<T>::operator!=(const QuaternionT<T>& q) const
{
return !(*this == q);
}

// Compute the quaternion that rotates from a to b, avoiding numerical instability.
// Taken from "The Shortest Arc Quaternion" by Stan Melax in "Game Programming Gems."
template <typename T>
inline QuaternionT<T> QuaternionT<T>::CreateFromVectors(const Vector3<T>& v0,
const Vector3<T>& v1)
{
if (v0 == -v1)
return QuaternionT<T>::CreateFromAxisAngle(vec3(1, 0, 0), Pi);

Vector3<T> c = v0.Cross(v1);
T d = v0.Dot(v1);
T s = std::sqrt((1 + d) * 2);

QuaternionT<T> q;
q.x = c.x / s;
q.y = c.y / s;
q.z = c.z / s;
q.w = s / 2.0f;
return q;
}

template <typename T>
inline QuaternionT<T>  QuaternionT<T>::CreateFromAxisAngle(const Vector3<T>& axis,
{
QuaternionT<T> q;
q.x = q.y = q.z = std::sin(radians / 2);
q.x *= axis.x;
q.y *= axis.y;
q.z *= axis.z;
return q;
}

template <typename T>
inline void QuaternionT<T>::Normalize()
{
*this = Scaled(1 / std::sqrt(Dot(*this)));
}

template <typename T>
inline void QuaternionT<T>::Rotate(const QuaternionT<T>& q2)
{
QuaternionT<T> q;
QuaternionT<T>& q1 = *this;

q.w = q1.w * q2.w - q1.x * q2.x - q1.y * q2.y - q1.z * q2.z;
q.x = q1.w * q2.x + q1.x * q2.w + q1.y * q2.z - q1.z * q2.y;
q.y = q1.w * q2.y + q1.y * q2.w + q1.z * q2.x - q1.x * q2.z;
q.z = q1.w * q2.z + q1.z * q2.w + q1.x * q2.y - q1.y * q2.x;

q.Normalize();
*this = q;
}

typedef QuaternionT<float> Quaternion;