I have a matrix class that allows me to perform rotation about an arbitrary axis. The problem is, I can't seem to figure out how to perform rotations relative to previous rotations. For example, if I rotate 90 degrees about the y-axis, I want the x-axis and z-axis to also rotate 90 degrees and affect future rotations done in this new coordinate space. How can I do this? In my naive attempt I tried to transform all the inputed rotation vectors by the matrix before rotating about the vector, this didn't work though.

Composition of rotations corresponds to matrix multiplication. If it doesn't seem to be behaving the way you expect, perhaps you are multiplying the matrices in the wrong order.

If that doesn't help you, please post a concrete example of two rotations that are not being compounded correctly.

Here is an example, it translates as expected, rotates about the y-axis as expected but then the x-axis rotation is done as if the y-axis was not rotated at all.

[source]
void doSomething() {
_matrix.loadIdentity();
_matrix.translate(Vector3f(100.0f, 200.0f, 0.0f));
_matrix.rotate(Vector3f(0.0f, 1.0f, 0.0f), 30.0f);
_matrix.rotate(Vector3f(1.0f, 0.0f, 0.0f), 45.0f);
}

class Matrix4x4f {
private:
float _data[16];

public:
void translate(const Vector3f& translation);
void rotate(const Vector3f& axis, float angle);
};

void Matrix4x4f::translate(const Vector3f& translation) {
float x = _data[0] * translation.x + _data[4] * translation.y + _data[8] * translation.z + _data[12];
float y = _data[1] * translation.x + _data[5] * translation.y + _data[9] * translation.z + _data[13];
float z = _data[2] * translation.x + _data[6] * translation.y + _data[10] * translation.z + _data[14];
float w = _data[3] * translation.x + _data[7] * translation.y + _data[11] * translation.z + _data[15];

_data[12] = x;
_data[13] = y;
_data[14] = z;
_data[15] = w;
}

void Matrix4x4f::rotate(const Vector3f& axis, float angle) {
float angleInRadians = BasicMath::toRadians(angle);
float cosAngle = cosf(angleInRadians);
float sinAngle = sinf(angleInRadians);
float cosDifference = 1.0f - cosAngle;

float rotationMatrix[9];
rotationMatrix[0] = cosDifference * axis.x * axis.x + cosAngle;
rotationMatrix[1] = cosDifference * axis.x * axis.y + sinAngle * axis.z;
rotationMatrix[2] = cosDifference * axis.x * axis.z - sinAngle * axis.y;

rotationMatrix[3] = cosDifference * axis.x * axis.y - sinAngle * axis.z;
rotationMatrix[4] = cosDifference * axis.y * axis.y + cosAngle;
rotationMatrix[5] = cosDifference * axis.y * axis.z + sinAngle * axis.x;

rotationMatrix[6] = cosDifference * axis.x * axis.z + sinAngle * axis.y;
rotationMatrix[7] = cosDifference * axis.y * axis.z - sinAngle * axis.x;
rotationMatrix[8] = cosDifference * axis.z * axis.z + cosAngle;

float originalMatrix[9];
originalMatrix[0] = _data[0];
originalMatrix[1] = _data[1];
originalMatrix[2] = _data[2];

originalMatrix[3] = _data[4];
originalMatrix[4] = _data[5];
originalMatrix[5] = _data[6];

originalMatrix[6] = _data[8];
originalMatrix[7] = _data[9];
originalMatrix[8] = _data[10];

_data[0] = originalMatrix[0] * rotationMatrix[0] + originalMatrix[3] * rotationMatrix[1] + originalMatrix[6] * rotationMatrix[2];
_data[1] = originalMatrix[1] * rotationMatrix[0] + originalMatrix[4] * rotationMatrix[1] + originalMatrix[7] * rotationMatrix[2];
_data[2] = originalMatrix[2] * rotationMatrix[0] + originalMatrix[5] * rotationMatrix[1] + originalMatrix[8] * rotationMatrix[2];

_data[4] = originalMatrix[0] * rotationMatrix[3] + originalMatrix[3] * rotationMatrix[4] + originalMatrix[6] * rotationMatrix[5];
_data[5] = originalMatrix[1] * rotationMatrix[3] + originalMatrix[4] * rotationMatrix[4] + originalMatrix[7] * rotationMatrix[5];
_data[6] = originalMatrix[2] * rotationMatrix[3] + originalMatrix[5] * rotationMatrix[4] + originalMatrix[8] * rotationMatrix[5];

_data[8] = originalMatrix[0] * rotationMatrix[6] + originalMatrix[3] * rotationMatrix[7] + originalMatrix[6] * rotationMatrix[8];
_data[9] = originalMatrix[1] * rotationMatrix[6] + originalMatrix[4] * rotationMatrix[7] + originalMatrix[7] * rotationMatrix[8];
_data[10] = originalMatrix[2] * rotationMatrix[6] + originalMatrix[5] * rotationMatrix[7] + originalMatrix[8] * rotationMatrix[8];
}
[/source]

I'd seriously recommend breaking up your functions such that matrix multiplication is independent. Then implement your other functions in terms of that matrix multiplication. Keep in mind that matrix multiplication is non-commutative and the order that they are composed makes a difference.

The code seems correct to me. The way you are composing the matrices, the operations will happen in the opposite order to what you specified.

I turned your code into something I could run (it would have been nice if you had done this):
#include <iostream>
#include <cmath>

struct Vector3f {
float x, y, z;
Vector3f(float x, float y, float z) : x(x), y(y), z(z) {
}
};

struct Matrix4x4f {
float data[16];

void translate(const Vector3f& translation);
void rotate(const Vector3f& axis, float angle);
void loadIdentity();
};

void Matrix4x4f::translate(const Vector3f& translation) {
float x = data[0] * translation.x + data[4] * translation.y + data[8] * translation.z + data[12];
float y = data[1] * translation.x + data[5] * translation.y + data[9] * translation.z + data[13];
float z = data[2] * translation.x + data[6] * translation.y + data[10] * translation.z + data[14];
float w = data[3] * translation.x + data[7] * translation.y + data[11] * translation.z + data[15];

data[12] = x;
data[13] = y;
data[14] = z;
data[15] = w;
}

static const float radians_per_degree = std::atan(1.0f)/45.0f;

void Matrix4x4f::rotate(const Vector3f& axis, float angle) {
float angleInRadians = angle*radians_per_degree;
float cosAngle = std::cos(angleInRadians);
float sinAngle = std::sin(angleInRadians);
float cosDifference = 1.0f - cosAngle;

float rotationMatrix[9];
rotationMatrix[0] = cosDifference * axis.x * axis.x + cosAngle;
rotationMatrix[1] = cosDifference * axis.x * axis.y + sinAngle * axis.z;
rotationMatrix[2] = cosDifference * axis.x * axis.z - sinAngle * axis.y;

rotationMatrix[3] = cosDifference * axis.x * axis.y - sinAngle * axis.z;
rotationMatrix[4] = cosDifference * axis.y * axis.y + cosAngle;
rotationMatrix[5] = cosDifference * axis.y * axis.z + sinAngle * axis.x;

rotationMatrix[6] = cosDifference * axis.x * axis.z + sinAngle * axis.y;
rotationMatrix[7] = cosDifference * axis.y * axis.z - sinAngle * axis.x;
rotationMatrix[8] = cosDifference * axis.z * axis.z + cosAngle;

float originalMatrix[9];
originalMatrix[0] = data[0];
originalMatrix[1] = data[1];
originalMatrix[2] = data[2];

originalMatrix[3] = data[4];
originalMatrix[4] = data[5];
originalMatrix[5] = data[6];

originalMatrix[6] = data[8];
originalMatrix[7] = data[9];
originalMatrix[8] = data[10];

data[0] = originalMatrix[0] * rotationMatrix[0] + originalMatrix[3] * rotationMatrix[1] + originalMatrix[6] * rotationMatrix[2];
data[1] = originalMatrix[1] * rotationMatrix[0] + originalMatrix[4] * rotationMatrix[1] + originalMatrix[7] * rotationMatrix[2];
data[2] = originalMatrix[2] * rotationMatrix[0] + originalMatrix[5] * rotationMatrix[1] + originalMatrix[8] * rotationMatrix[2];

data[4] = originalMatrix[0] * rotationMatrix[3] + originalMatrix[3] * rotationMatrix[4] + originalMatrix[6] * rotationMatrix[5];
data[5] = originalMatrix[1] * rotationMatrix[3] + originalMatrix[4] * rotationMatrix[4] + originalMatrix[7] * rotationMatrix[5];
data[6] = originalMatrix[2] * rotationMatrix[3] + originalMatrix[5] * rotationMatrix[4] + originalMatrix[8] * rotationMatrix[5];

data[8] = originalMatrix[0] * rotationMatrix[6] + originalMatrix[3] * rotationMatrix[7] + originalMatrix[6] * rotationMatrix[8];
data[9] = originalMatrix[1] * rotationMatrix[6] + originalMatrix[4] * rotationMatrix[7] + originalMatrix[7] * rotationMatrix[8];
data[10] = originalMatrix[2] * rotationMatrix[6] + originalMatrix[5] * rotationMatrix[7] + originalMatrix[8] * rotationMatrix[8];
}

void Matrix4x4f::loadIdentity() {
for (int i=0; i<4; ++i) {
for (int j=0; j<4; ++j)
data[4*i+j] = (i==j);
}
}

Vector3f apply_to_point(Matrix4x4f const &m, Vector3f const &p) {
float v[4];
v[0] = p.x;
v[1] = p.y;
v[2] = p.z;
v[3] = 1.0f;
float r[4];
for (int i=0; i<4; ++i) {
float sum = 0.0f;
for (int j=0; j<4; ++j)
sum += m.data[i+4*j]*v[j];
r = sum;
}
return Vector3f(r[0],r[1],r[2]);
}

// Maps the point (1,2,3) through matrix m and prints the result
void display(Matrix4x4f const &m) {
Vector3f p(1.0f, 2.0f, 3.0f);
Vector3f r = apply_to_point(m,p);
std::cout << '(' << r.x << ',' << r.y << ',' << r.z << ")\n\n";
}

int main() {
Matrix4x4f matrix;
matrix.loadIdentity();
display(matrix);

matrix.translate(Vector3f(100.0f, 200.0f, 0.0f));
display(matrix);

matrix.rotate(Vector3f(0.0f, 1.0f, 0.0f), 30.0f);
display(matrix);

matrix.rotate(Vector3f(1.0f, 0.0f, 0.0f), 45.0f);
display(matrix);
}



Thanks for the advice guys and thanks for taking the time to write all that code Alvaro... sorry I didn't post a runnable sample. Anyway, it turns out that it's definitely the order of my transformations that was screwing me up. I'll have to fix my multiplication. Thanks again.