Simplifying Matrix Math

• Home
• Articles
• Contents
• Overview
• A Common Libary
• Simplifying Matrix Math
• The Code
• Performance Results

float f[16] =
{
    1.0f, 0.0f, 0.0f, x,
    0.0f, 1.0f, 0.0f, y,
    0.0f, 0.0f, 1.0f, z,
    0.0f, 0.0f, 0.0f, 1.0f
};

float f[16] =
{
    1.0f, 0.0f, 0.0f, 0.0f,
    0.0f, 1.0f, 0.0f, 0.0f,
    0.0f, 0.0f, 1.0f, 0.0f,
    x,    y,    z,    1.0f
};

I'm not going to work out all of the matrix simplification here. It would be much too cumbersome to format it in HTML. But do have a look at the library of code in this zip file: crossplat.zip. This is the same library for the Windows/MacOS compatibility articles I have been writing. The relevant file is QzMatrix4x4.h and cpp.

For a quick example, I'll show a simple translation operation. Let's assume that you start with a transform matrix filled with unknown coefficients, each represented by the letters a–p.

float matrix[16] =
{
    a, b, c, d,
    e, f, g, h,
    i, j, k, l,
    m, n, o, p
};

If you want to apply a translation by (x,y,z), you would construct this matrix.

float translate[16] =
{
    1, 0, 0, 0,
    0, 1, 0, 0,
    0, 0, 1, 0,
    x, y, z, 1
};

To apply the translation such that it occurs before any other transformation, you would perform the matrix multiplication translate × matrix (unless, of course, you've been taught to multiply matrices the other way around, in which case it would be matrix × translate).

The result of this 4×4 matrix multiplication would be the following matrix:

float f[16] =
{
    a, b, c, d,
    e, f, g, h,
    i, j, k, l,
    x*a + y*e + z*i + m,
    x*b + y*f + z*j + n,
    x*c + y*g + z*k + o,
    x*d + y*h + z*l + p
};

Note that only the bottom row of coefficients is changed. The other values are completely left alone. Similarly, all other variations of translation, scaling, and rotation can be simplified, leaving many coefficients unchanged.

By factoring the full 4×4 matrix multiplication out of the picture, all of the standard transformations will run much faster. The only expensive operation that remains is for rotations, which must still compute a sin and cos — indeed, these two transcendental functions combined can be much more expensive to execute than a full 4×4 matrix multiplication, depending upon the CPU type (see the results at the end of the article).

Many of the routines in my library are prefixed by Pre... or Post.... Pre-transforms are applied before any other operation, and post-transforms are applied afterwards. For example, if you needed to scale a model down 50%, then rotate it around the X-axis, and finally translate it to a specific (x,y,z) position, you could do it with the following sequence of operations:

    matrix.MakeScale(0.5f);
    matrix.PostRotateX(angle);
    matrix.PostTranslate(x, y, z);

On the other hand, you may be accustomed to building up a transform matrix in the reverse order, starting with the last transform and working down to to first transform:

    matrix.MakeTranslate(x, y, z);
    matrix.PreRotateX(angle);
    matrix.PreScale(0.5f);

I've found that explicitly naming the transforms as pre- or post- eliminates the confusion surrounding matrix composition. Mileage may vary. Your very notion of "pre" and "post" may be the exact opposite of mine.

struct QzMatrix4x4_t
{
    union {
        float c4x4[4][4];
        float cFlat[16];
    };
};

For anyone not having access to original source file, QzMatrix4x4.cpp, here is a copy of the relevant code:

void QzMatrix4x4::Multiply(const QzMatrix4x4 &first, const QzMatrix4x4 &second)
{
    const float *m1 = first.m_Matrix.cFlat;
    const float *m2 = second.m_Matrix.cFlat;
          float *r  = m_Matrix.cFlat;

    for (U32 row = 0; row < 4; ++row) {
        for (U32 col = 0; col < 4; ++col) {
            r[(row<<2)+col] = m1[(row<<2)  ] * m2[col   ]
                            + m1[(row<<2)+1] * m2[col+ 4]
                            + m1[(row<<2)+2] * m2[col+ 8]
                            + m1[(row<<2)+3] * m2[col+12];
        }
    }
}

void QzMatrix4x4::PreMultiply(const QzMatrix4x4 &first)
{
    // Multiplication is destructive, so we need make a duplicate of the
    // original value of m_Matrix for use in computing the new m_Matrix.
    const QzMatrix4x4_t dup = m_Matrix;
    const float*        m1  = first.m_Matrix.cFlat;
    const float*        m2  = dup.cFlat;
          float*        r   = m_Matrix.cFlat;

    for (U32 row = 0; row < 4; ++row) {
        for (U32 col = 0; col < 4; ++col) {
            r[(row<<2)+col] = m1[(row<<2)  ] * m2[col   ]
                            + m1[(row<<2)+1] * m2[col+ 4]
                            + m1[(row<<2)+2] * m2[col+ 8]
                            + m1[(row<<2)+3] * m2[col+12];
        }
    }
}

void QzMatrix4x4::PostMultiply(const QzMatrix4x4 &second)
{
    // Multiplication is destructive, so we need make a duplicate of the
    // original value of m_Matrix for use in computing the new m_Matrix.
    const QzMatrix4x4_t dup = m_Matrix;
    const float*        m1  = dup.cFlat;
    const float*        m2  = second.m_Matrix.cFlat;
          float*        r   = m_Matrix.cFlat;

    for (U32 row = 0; row < 4; ++row) {
        for (U32 col = 0; col < 4; ++col) {
            r[(row<<2)+col] = m1[(row<<2)  ] * m2[col   ]
                            + m1[(row<<2)+1] * m2[col+ 4]
                            + m1[(row<<2)+2] * m2[col+ 8]
                            + m1[(row<<2)+3] * m2[col+12];
        }
    }
}

void QzMatrix4x4::PreRotateX(float angle)
{
    float c = cosf(angle);
    float s = sinf(angle);

    float *p = m_Matrix.cFlat;

    for (U32 i = 0; i < 4; ++i) {
        float t = p[4+i];
        float u = p[8+i];
        p[4+i]  = (t * c) + (u * s);
        p[8+i]  = (u * c) - (t * s);
    }
}

void QzMatrix4x4::PostRotateX(float angle)
{
    float c = cosf(angle);
    float s = sinf(angle);

    float *p = m_Matrix.cFlat;

    for (U32 i = 0; i < 4; ++i) {
        float t = p[(i<<2)+1];
        float u = p[(i<<2)+2];
        p[(i<<2)+1] = (t * c) - (u * s);
        p[(i<<2)+2] = (t * s) + (u * c);
    }
}

void QzMatrix4x4::PreRotateY(float angle)
{
    float c = cosf(angle);
    float s = sinf(angle);

    float *p = m_Matrix.cFlat;

    for (U32 i = 0; i < 4; ++i) {
        float t = p[  i];
        float u = p[8+i];
        p[  i]  = (t * c) - (u * s);
        p[8+i]  = (u * c) + (t * s);
    }
}

void QzMatrix4x4::PostRotateY(float angle)
{
    float c = cosf(angle);
    float s = sinf(angle);

    float *p = m_Matrix.cFlat;

    for (U32 i = 0; i < 4; ++i) {
        float t = p[(i<<2)  ];
        float u = p[(i<<2)+2];
        p[(i<<2)  ] = (u * s) + (t * c);
        p[(i<<2)+2] = (u * c) - (t * s);
    }
}

void QzMatrix4x4::PreRotateZ(float angle)
{
    float c = cosf(angle);
    float s = sinf(angle);

    float *p = m_Matrix.cFlat;

    for (U32 i = 0; i < 4; ++i) {
        float t = p[  i];
        float u = p[4+i];
        p[  i]  = (t * c) + (u * s);
        p[4+i]  = (u * c) - (t * s);
    }
}

void QzMatrix4x4::PostRotateZ(float angle)
{
    float c = cosf(angle);
    float s = sinf(angle);

    float *p = m_Matrix.cFlat;

    for (U32 i = 0; i < 4; ++i) {
        float t = p[(i<<2)  ];
        float u = p[(i<<2)+1];
        p[(i<<2)  ] = (t * c) - (u * s);
        p[(i<<2)+1] = (t * s) + (u * c);
    }
}

void QzMatrix4x4::PreScale(float x, float y, float z)
{
    for (U32 i = 0; i < 4; ++i) {
        m_Matrix.c4x4[0][i] *= x;
        m_Matrix.c4x4[1][i] *= y;
        m_Matrix.c4x4[2][i] *= z;
    }
}

void QzMatrix4x4::PreScale(float x)
{
    for (U32 i = 0; i < 4; ++i) {
        m_Matrix.c4x4[0][i] *= x;
        m_Matrix.c4x4[1][i] *= x;
        m_Matrix.c4x4[2][i] *= x;
    }
}

void QzMatrix4x4::PostScale(float x, float y, float z)
{
    for (U32 i = 0; i < 4; ++i) {
        m_Matrix.c4x4[i][0] *= x;
        m_Matrix.c4x4[i][1] *= y;
        m_Matrix.c4x4[i][2] *= z;
    }
}

void QzMatrix4x4::PostScale(float x)
{
    for (U32 i = 0; i < 4; ++i) {
        m_Matrix.c4x4[i][0] *= x;
        m_Matrix.c4x4[i][1] *= x;
        m_Matrix.c4x4[i][2] *= x;
    }
}

void QzMatrix4x4::PreTranslate(float x, float y, float z)
{
    for (U32 i = 0; i < 4; ++i) {
        m_Matrix.c4x4[3][i] += (x * m_Matrix.c4x4[0][i])
                            +  (y * m_Matrix.c4x4[1][i])
                            +  (z * m_Matrix.c4x4[2][i]);
    }
}

void QzMatrix4x4::PostTranslate(float x, float y, float z)
{
    for (U32 i = 0; i < 4; ++i) {
        m_Matrix.c4x4[i][0] += m_Matrix.c4x4[i][3] * x;
        m_Matrix.c4x4[i][1] += m_Matrix.c4x4[i][3] * y;
        m_Matrix.c4x4[i][2] += m_Matrix.c4x4[i][3] * z;
    }
}

The zip file, crossplat.zip, contains a QzUnitTest project. Running this project will verify that all of these transforms yield the same results as the traditional matrix multiplication approach.

The unit test in the zip file will also perform some benchmarking of these operations. The results are written out to trace.txt. I've benchmarked a few machines to give an idea how well these operations run.

On Windows, these numbers were taken from code compiled with the 2003 Professional version of DevStudio, which has an optimizing compiler. The Standard and Express versions of DevStudio seldom (never?) have optimizing compilers, so tightly written C code like matrix multiplication can easily take several times as many CPU cycles to execute as code generated by an optimizing compiler.

As an example of how much difference this makes, consider the Multiply function. On a version of DevStudio with an optimizing compiler, a release build takes 0.057 microseconds, or 0.299 microseconds with a debug build. Compiling the same code on the Standard 2003 version of DevStudio, this function takes 0.212 microseconds, almost four times as long as code generated by an optimizing compiler. Obviously, if you care about performance on Windows, you will have to spring the several hundred dollars required for a Professional version of DevStudio.

So what do the numbers in the previous table indicate?

Going by the numbers from the Core 2 6600, translating an object would require a call to MakeTranslate to create the 4×4 translation matrix, then a call to MatrixMultiply, for a total of 0.057 + 0.016 = 0.073 μs. Directly applying the transform via PostTranslate would only take 0.017 μs, or about a 4× difference in processing time.

Creating a rotation matrix would require 0.083 μs, or 0.140 μs counting the matrix multiply, versus the 0.091 μs required to apply the rotation directly with PostRotateX. Here the difference is not as pronounced, since there is a 0.073 μs overhead required for the sinf and cosf function calls needed to compute the rotation coefficients. Still, most of the sampled processors exhibit a 30% to 50% reduction in CPU time for rotations. And when running with optimized sine functions, the performance improvement is even larger.

In some cases, the pre- and post- transforms require copying the contents of the matrix to avoid order dependence problems with computing new coefficients. It may only involve copying 16 floats, but this will still add a few nanoseconds to the transform time.

How well do these simplifications perform on embedded processors? That heavily depends on CPU speed versus memory speed. When the CPU is significantly faster than memory, a tight matrix multiply loop gains a preformance improvement since it requires a smaller instruction cache footprint — the simplified transforms should still be faster, but the difference may not be very dramatic. On the other hand, when CPU speed is closer to memory speed (often the case with embedded devices), there is a greater performance improvement when the total number of CPU instructions is reduced, so the simplifications in this article provide a greater benefit — even though the memory footprint of the unrolled instructions is larger than a tight loop.

With today's PC hardware, the speed of matrix operations is not likely to be a cost factor (unless you're doing something like skinned animations in software), but the weak embedded processors used in handhelds can still benefit from every clock cycle that is saved by simplifying the matrix transforms.