4D matrix math library for use with OpenGL

Question

I am trying to create a simple library for C to handle OpenGL matrix operations.

You can see the vec3fscalar here.

#ifndef __MATH_MAT4FSCALAR_H__
#define __MATH_MAT4FSCALAR_H__

#include "vec3fScalar.h"

typedef struct mat4
{
    double m[16];
}mat4;

mat4* newMat4(const vec3f *c1, const vec3f *c2, const vec3f *c3, const vec3f *c4);
mat4* newMat4identity(void);

void newMat4x(mat4* out, const vec3f *c1, const vec3f *c2, const vec3f *c3, const vec3f *c4);
void newMat4xc(mat4* out, const mat4* m4);
void newMat4identityx(mat4* out);


void mat4add(mat4* out, const mat4* a, const mat4* b);
void mat4sub(mat4* out, const mat4* a, const mat4* b);
void mat4mul(mat4* out, const mat4* a, const mat4* b);

void mat4mulvec3(vec3f* out, const vec3f* v, const mat4* a);

void mat4trans(mat4* out, const double x, const double y, const double z);
void mat4transv(mat4* out, const vec3f* v);

void mat4rotate(mat4* out, const double angle, const double x, const double y, const double z);
void mat4rotatev(mat4* out, const double angle, const vec3f* v);
void mat4xrotate(mat4* out, const double angle, const double x);
void mat4yrotate(mat4* out, const double angle, const double y);
void mat4zrotate(mat4* out, const double angle, const double z);


void mat4inv(mat4* out, const mat4* in);

void mat4print(const char* tag, const mat4* m4);

#endif /* __MATH_AT4FSCALAR_H__ */

mat4.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "mat4fScalar.h"

static double iData[] = { 1, 0, 0, 0,
                          0, 1, 0, 0,
                          0, 0, 1, 0,
                          0, 0, 0, 1 };

static mat4* mMat(mat4* out, const double* data)
{
    out->m[0]  = data[0];
    out->m[4]  = data[4];
    out->m[8]  = data[8];
    out->m[12] = data[12];

    out->m[1]  = data[1];
    out->m[5]  = data[5];
    out->m[9]  = data[9];
    out->m[13] = data[13];

    out->m[2]  = data[2];
    out->m[6]  = data[6];
    out->m[10] = data[10];
    out->m[14] = data[14];

    out->m[3]  = data[3];
    out->m[7]  = data[7];
    out->m[11] = data[11];
    out->m[15] = data[15];
    return out;

}

mat4* newMat4(const vec3f *c1, const vec3f *c2, const vec3f *c3, const vec3f *c4)
{
    double data[] = { c1->p[0], c2->p[0], c3->p[0], c4->p[0],
                      c1->p[1], c2->p[1], c3->p[1], c4->p[1],
                      c1->p[2], c2->p[2], c3->p[2], c4->p[2],
                         0,        0,       0,        1      };

    mat4* v = (mat4*)malloc(sizeof(mat4));
    return mMat(v, data);
}

mat4* newMat4c(const mat4* m4)
{
    mat4* v = (mat4*)malloc(sizeof(mat4));
    return mMat(v, m4->m);
}

mat4* newMat4identity(void)
{
    mat4* v = (mat4*)malloc(sizeof(mat4));
    return mMat(v, iData);
}

void newMat4x(mat4* out, const vec3f *c1, const vec3f *c2, const vec3f *c3, const vec3f *c4)
{
    double data[] = { c1->p[0], c2->p[0], c3->p[0], c4->p[0],
                      c1->p[1], c2->p[1], c3->p[1], c4->p[1],
                      c1->p[2], c2->p[2], c3->p[2], c4->p[2],
                      0,        0,       0,        1      };
     mMat(out, data);
}

void newMat4xc(mat4* out, const mat4* m4)
{
    mMat(out, m4->m);
}

void newMat4identityx(mat4* out)
{
    mMat(out, iData);
}

void mat4add(mat4* out, const mat4* a, const mat4* b)
{
    double d[16];
    for(size_t i = 0; i < 16; i++)
    {
        d[i] = a->m[i] + b->m[i];
    }
    mMat(out, d);
}

void mat4sub(mat4* out, const mat4* a, const mat4* b)
{
    double d[16];
    for(size_t i = 0; i < 16; i++)
    {
        d[i] = a->m[i] - b->m[i];
    }
    mMat(out, d);
}

void mat4mul(mat4* out, const mat4* a, const mat4* b)
{
    out->m[0]  = ( a->m[0]  *  b->m[0]) + (a->m[1]  * b->m[4]) + (a->m[2]  *  b->m[8]) + (a->m[3]  * b->m[12] );
    out->m[4]  = ( a->m[4]  *  b->m[0]) + (a->m[5]  * b->m[4]) + (a->m[6]  *  b->m[8]) + (a->m[7]  * b->m[12] );
    out->m[8]  = ( a->m[8]  *  b->m[0]) + (a->m[9]  * b->m[4]) + (a->m[10] *  b->m[8]) + (a->m[11] * b->m[12] );
    out->m[12] = ( a->m[12] *  b->m[0]) + (a->m[13] * b->m[4]) + (a->m[14] *  b->m[8]) + (a->m[15] * b->m[12] );

    out->m[1]  = ( a->m[0]  *  b->m[1]) + (a->m[1]  * b->m[5]) + (a->m[2]  * b->m[9])  + (a->m[3]  * b->m[13] );
    out->m[5]  = ( a->m[4]  *  b->m[1]) + (a->m[5]  * b->m[5]) + (a->m[6]  * b->m[9])  + (a->m[7]  * b->m[13] );
    out->m[9]  = ( a->m[8]  *  b->m[1]) + (a->m[9]  * b->m[5]) + (a->m[10] * b->m[9])  + (a->m[11] * b->m[13] );
    out->m[13] = ( a->m[12] *  b->m[1]) + (a->m[13] * b->m[5]) + (a->m[14] * b->m[9])  + (a->m[15] * b->m[13] );

    out->m[2]  = ( a->m[0]  *  b->m[2]) + (a->m[1]  * b->m[6]) + (a->m[2]  * b->m[10]) + (a->m[3]  * b->m[14] );
    out->m[6]  = ( a->m[4]  *  b->m[2]) + (a->m[5]  * b->m[6]) + (a->m[6]  * b->m[10]) + (a->m[7]  * b->m[14] );
    out->m[10] = ( a->m[8]  *  b->m[2]) + (a->m[9]  * b->m[6]) + (a->m[10] * b->m[10]) + (a->m[11] * b->m[14] );
    out->m[14] = ( a->m[12] *  b->m[2]) + (a->m[13] * b->m[6]) + (a->m[14] * b->m[10]) + (a->m[15] * b->m[14] );

    out->m[3]  = ( a->m[0]  *  b->m[3]) + (a->m[1]  * b->m[7]) + (a->m[2]  * b->m[11]) + (a->m[3]  * b->m[15] );
    out->m[7]  = ( a->m[4]  *  b->m[3]) + (a->m[5]  * b->m[7]) + (a->m[6]  * b->m[11]) + (a->m[7]  * b->m[15] );
    out->m[11] = ( a->m[8]  *  b->m[3]) + (a->m[9]  * b->m[7]) + (a->m[10] * b->m[11]) + (a->m[11] * b->m[15] );
    out->m[15] = ( a->m[12] *  b->m[3]) + (a->m[13] * b->m[7]) + (a->m[14] * b->m[11]) + (a->m[15] * b->m[15] );
}

void mat4mulvec3(vec3f* out, const vec3f* v, const mat4* a)
{
    out->p[0] = a->m[0] * v->p[0] + a->m[1] * v->p[1] + a->m[2]  * v->p[2];
    out->p[1] = a->m[4] * v->p[0] + a->m[5] * v->p[1] + a->m[6]  * v->p[2];
    out->p[2] = a->m[8] * v->p[0] + a->m[9] * v->p[1] + a->m[10] * v->p[2];
}

void mat4trans(mat4* out, const double x, const double y, const double z)
{
    out->m[3]  += x;
    out->m[7]  += y;
    out->m[11] += z;
}

void mat4transv(mat4* out, const vec3f* v)
{
    out->m[3]  += v->p[0];
    out->m[7]  += v->p[1];
    out->m[11] += v->p[2];
}

void mat4rotate(mat4*  out, const double angle, const double x, const double y, const double z)
{
    double mAng = toRadiansr(angle);
    out->m[0]   = cos(y) * cos(z * mAng);
    out->m[4]   = cos(y) * sin(z * mAng);
    out->m[8]   = sin(y);

    out->m[1]   = ((cos(y) * sin(z)) + (sin(x) * sin(y) * cos(z)));
    out->m[5]   = ((cos(y) * cos(z)) - (sin(x) * sin(y) * sin(z)));
    out->m[9]   = -sin(x)  * cos(y);

    out->m[2]   = ((sin(x) * sin(z)) - (cos(x) * sin(y) * cos(z)));
    out->m[6]   = ((sin(x) * sin(z)) + (cos(x) * sin(y) * sin(z)));
    out->m[10]  = (cos(x)  * cos(y));
}

void mat4rotatev(mat4* out, const double angle, const vec3f* v)
{
    double mAng = toRadiansr(angle);
    double x, y, z;
    x = v->p[0];
    y = v->p[1];
    z = v->p[2];

    out->m[0]   = cos(y) * cos(z * mAng);
    out->m[4]   = cos(y) * sin(z * mAng);
    out->m[8]   = sin(y);

    out->m[1]   = ((cos(y) * sin(z)) + (sin(x) * sin(y) * cos(z)));
    out->m[5]   = ((cos(y) * cos(z)) - (sin(x) * sin(y) * sin(z)));
    out->m[9]   = -sin(x)  * cos(y);

    out->m[2]   = ((sin(x) * sin(z)) - (cos(x) * sin(y) * cos(z)));
    out->m[6]   = ((sin(x) * sin(z)) + (cos(x) * sin(y) * sin(z)));
    out->m[10]  = (cos(x)  * cos(y));
}


void mat4xrotate(mat4* out, const double angle, const double x)
{
    double mAng = toRadiansr(angle);
    out->m[5]  = cos(mAng);
    out->m[9]  = -sin(mAng);

    out->m[6]  = sin(mAng);
    out->m[10] = cos(mAng);
}

void mat4yrotate(mat4* out, const double angle, const double y)
{
    double mAng = toRadiansr(angle);
    out->m[0]  = cos(mAng);
    out->m[8]  = sin(mAng);

    out->m[2]  = -sin(mAng);
    out->m[10] = cos(mAng);
}

void mat4zrotate(mat4* out, const double angle, const double z)
{
    double mAng = toRadiansr(angle);
    out->m[0] = cos(mAng);
    out->m[4] = sin(mAng);

    out->m[1] = -sin(mAng);
    out->m[5] = cos(mAng);
}

void mat4scale(mat4* out, const double x, const double y, const double z)
{
    out->m[0]   *= x;
    out->m[5]   *= y;
    out->m[10]  *= z;
}

void mat4scalev(mat4* out, const vec3f* v)
{
    out->m[0]   *= v->p[0];
    out->m[5]   *= v->p[1];
    out->m[10]  *= v->p[2];
}

void mat4inv(mat4* out, const mat4* in)
{
    double inv[16], det;
    int i;

    inv[0] =  in->m[5]  * in->m[10] * in->m[15] -
              in->m[5]  * in->m[11] * in->m[14] -
              in->m[9]  * in->m[6]  * in->m[15] +
              in->m[9]  * in->m[7]  * in->m[14] +
              in->m[13] * in->m[6]  * in->m[11] -
              in->m[13] * in->m[7]  * in->m[10];

    inv[4] = -in->m[4]  * in->m[10] * in->m[15] +
              in->m[4]  * in->m[11] * in->m[14] +
              in->m[8]  * in->m[6]  * in->m[15] -
              in->m[8]  * in->m[7]  * in->m[14] -
              in->m[12] * in->m[6]  * in->m[11] +
              in->m[12] * in->m[7]  * in->m[10];

    inv[8] =  in->m[4]  * in->m[9]  * in->m[15] -
              in->m[4]  * in->m[11] * in->m[13] -
              in->m[8]  * in->m[5]  * in->m[15] +
              in->m[8]  * in->m[7]  * in->m[13] +
              in->m[12] * in->m[5]  * in->m[11] -
              in->m[12] * in->m[7]  * in->m[9];

    inv[12] = -in->m[4]  * in->m[9] * in->m[14] +
               in->m[4]  * in->m[10] * in->m[13] +
               in->m[8]  * in->m[5] * in->m[14] -
               in->m[8]  * in->m[6] * in->m[13] -
               in->m[12] * in->m[5] * in->m[10] +
               in->m[12] * in->m[6] * in->m[9];

    inv[1] =  -in->m[1]  * in->m[10] * in->m[15] +
              in->m[1]  * in->m[11] * in->m[14] +
              in->m[9]  * in->m[2] * in->m[15] -
              in->m[9]  * in->m[3] * in->m[14] -
              in->m[13] * in->m[2] * in->m[11] +
              in->m[13] * in->m[3] * in->m[10];

    inv[5] =  in->m[0]  * in->m[10] * in->m[15] -
              in->m[0]  * in->m[11] * in->m[14] -
              in->m[8]  * in->m[2] * in->m[15] +
              in->m[8]  * in->m[3] * in->m[14] +
              in->m[12] * in->m[2] * in->m[11] -
              in->m[12] * in->m[3] * in->m[10];

    inv[9] =  -in->m[0]  * in->m[9] * in->m[15] +
              in->m[0]  * in->m[11] * in->m[13] +
              in->m[8]  * in->m[1] * in->m[15] -
              in->m[8]  * in->m[3] * in->m[13] -
              in->m[12] * in->m[1] * in->m[11] +
              in->m[12] * in->m[3] * in->m[9];

    inv[13] =  in->m[0]  * in->m[9] * in->m[14] -
               in->m[0]  * in->m[10] * in->m[13] -
               in->m[8]  * in->m[1] * in->m[14] +
               in->m[8]  * in->m[2] * in->m[13] +
               in->m[12] * in->m[1] * in->m[10] -
               in->m[12] * in->m[2] * in->m[9];

    inv[2] =  in->m[1]  * in->m[6] * in->m[15] -
              in->m[1]  * in->m[7] * in->m[14] -
              in->m[5]  * in->m[2] * in->m[15] +
              in->m[5]  * in->m[3] * in->m[14] +
              in->m[13] * in->m[2] * in->m[7] -
              in->m[13] * in->m[3] * in->m[6];

    inv[6] =  -in->m[0]  * in->m[6] * in->m[15] +
              in->m[0]  * in->m[7] * in->m[14] +
              in->m[4]  * in->m[2] * in->m[15] -
              in->m[4]  * in->m[3] * in->m[14] -
              in->m[12] * in->m[2] * in->m[7] +
              in->m[12] * in->m[3] * in->m[6];

    inv[10] =  in->m[0]  * in->m[5] * in->m[15] -
               in->m[0]  * in->m[7] * in->m[13] -
               in->m[4]  * in->m[1] * in->m[15] +
               in->m[4]  * in->m[3] * in->m[13] +
               in->m[12] * in->m[1] * in->m[7] -
               in->m[12] * in->m[3] * in->m[5];

    inv[14] =  -in->m[0]  * in->m[5] * in->m[14] +
               in->m[0]  * in->m[6] * in->m[13] +
               in->m[4]  * in->m[1] * in->m[14] -
               in->m[4]  * in->m[2] * in->m[13] -
               in->m[12] * in->m[1] * in->m[6] +
               in->m[12] * in->m[2] * in->m[5];

    inv[3] =  -in->m[1] * in->m[6] * in->m[11] +
              in->m[1] * in->m[7] * in->m[10] +
              in->m[5] * in->m[2] * in->m[11] -
              in->m[5] * in->m[3] * in->m[10] -
              in->m[9] * in->m[2] * in->m[7] +
              in->m[9] * in->m[3] * in->m[6];

    inv[7] =  in->m[0] * in->m[6] * in->m[11] -
              in->m[0] * in->m[7] * in->m[10] -
              in->m[4] * in->m[2] * in->m[11] +
              in->m[4] * in->m[3] * in->m[10] +
              in->m[8] * in->m[2] * in->m[7] -
              in->m[8] * in->m[3] * in->m[6];

    inv[11] = -in->m[0] * in->m[5] * in->m[11] +
              in->m[0] * in->m[7] * in->m[9] +
              in->m[4] * in->m[1] * in->m[11] -
              in->m[4] * in->m[3] * in->m[9] -
              in->m[8] * in->m[1] * in->m[7] +
              in->m[8] * in->m[3] * in->m[5];

    inv[15] = in->m[0] * in->m[5] * in->m[10] -
              in->m[0] * in->m[6] * in->m[9] -
              in->m[4] * in->m[1] * in->m[10] +
              in->m[4] * in->m[2] * in->m[9] +
              in->m[8] * in->m[1] * in->m[6] -
              in->m[8] * in->m[2] * in->m[5];

    det = in->m[0] * inv[0] + in->m[1] * inv[4] + in->m[2] * inv[8] + in->m[3] * inv[12];

    if (det == 0)
        printf("oh oh!");

    det = 1.0 / det;

    for (i = 0; i < 16; i++)
        out->m[i] = inv[i] * det;
}

void mat4print(const char* tag, const mat4* m4)
{
            printf("\n\t%s\t\n[\t%f\t%f\t%f\t%f\t]\n[\t%f\t%f\t%f\t%f\t]\n"
                           "[\t%f\t%f\t%f\t%f\t]\n[\t%f\t%f\t%f\t%f\t]\n",
                   tag,
                   m4->m[0],  m4->m[1],  m4->m[2],  m4->m[3],
                   m4->m[4],  m4->m[5],  m4->m[6],  m4->m[7],
                   m4->m[8],  m4->m[9],  m4->m[10], m4->m[11],
                   m4->m[12], m4->m[13], m4->m[14], m4->m[15]);
}

I didn't write my own inverse function; I got it from this SO question.

I really hope that I can get some good help with this because matrix math is difficult enough w/o translating it to code. Are there any functions that I am missing? I would like to avoid adding functions just to have them there unless they'll be useful for rendering pipeline which fits nicely in the 4x4 square matrix.

Also, please no comments about rolling up or using loops.

Answer 1

A few things you might want to look into:

• Prefer to avoid names starting with double underscore (__MATH_MAT4FSCALAR_H__). This notation is usually reserved for compiler extensions and Standard Library internals. Steer clear from them to avoid a name collision with your own code. In this case, simple MATH_MAT4FSCALAR_H would be just fine.

• I don't see any reason for not making iData a mat4 as well. Then also take two matrices in mMat(). By the way, mMat tells me nothing. copyMatrix would be the logical name for that function if it took two matrices and copied b over a. If you leave the second param as double array, then perhaps setMatrix will make more sense.

• Avoid casting the return value of malloc, unless you need to compile the code as C++. In C the cast from void* to whatever is implicit, so that only adds boilerplate and maintenance overhead.

• You might also consider using the variable instance in the sizeof on mallocs to avoid repeating the types on the left and right, E.g.:

  mat4* v = malloc(sizeof *v); // The ( ) around sizeof are not even mandatory!
  

• Descent compilers are capable of consolidating those redundant sin/cos calls on the same value during the common subexpression elimination pass, but... I don't know, call me paranoid, but I'd rather not rely on that. Specially if you target some more exotic systems with much simpler compilers. I'd manually consolidate the repeated calls into local variables and use the locals instead.

• Hehe, okay, suppose that does happen:

  if (det == 0)
      printf("oh oh!");
  

  What is your course of action? Allow the division by zero to take place, returning NaN/Inf? Better would be to make the function return a boolean and let the caller know about the error.

• In mat4print(), you might consider accepting a third parameter with the FILE* to print to. That makes you code reusable if you want for instance write to a log file on disk. For console dumps, just pass stdout.

---

Apart from the above, it looks good to me. I commend you for making consistent use of const. I'm a "const nazi" myself. Good job!