felselva / mathc

Pure C math library for 2D and 3D programming
zlib License
711 stars 54 forks source link

Adding inverse and transpose functions for the matrix. #11

Closed Ankush-p closed 6 years ago

Ankush-p commented 6 years ago

I was in need of two functions, one to transpose a matrix (trun it's columns into rows) and another function to invert the matrix. I needed these functions for OpenGL - after writing up these two functions, I decided to open the issue and place the code below, incase you (@ferreiradaselva) wish to add it to the library.

Transpose Matrix (pmatrix_transpose() and matrix_transpose())

void pmatrix_transpose(struct matrix *m, struct matrix *result)
{
    result->m11 = m->m11;
    result->m21 = m->m12;
    result->m31 = m->m13;
    result->m41 = m->m14;

    result->m12 = m->m21;
    result->m22 = m->m22;
    result->m32 = m->m23;
    result->m42 = m->m24;

    result->m13 = m->m31;
    result->m23 = m->m32;
    result->m33 = m->m33;
    result->m43 = m->m34;

    result->m14 = m->m41;
    result->m24 = m->m42;
    result->m34 = m->m43;
    result->m44 = m->m44;
}

struct matrix matrix_transpose(struct matrix m)
{
    struct matrix result;
    pmatrix_transpose(&m, &result);
    return (result);
}

Inverse Matrix (pmatrix_inverse() and matrix_inverse())

void pmatrix_inverse(struct matrix *m, struct matrix *result)
{
    struct matrix inv;
    float det;

    inv.m11 = m->m22 * m->m33 * m->m44 -
          m->m22 * m->m43 * m->m34 -
          m->m23 * m->m32 * m->m44 +
          m->m23 * m->m42 * m->m34 +
          m->m24 * m->m32 * m->m43 -
          m->m24 * m->m42 * m->m33;

    inv.m12 = -m->m12 * m->m33 * m->m44 +
          m->m12 * m->m43 * m->m34 +
          m->m13 * m->m32 * m->m44 -
          m->m13 * m->m42 * m->m34 -
          m->m14 * m->m32 * m->m43 +
          m->m14 * m->m42 * m->m33;

    inv.m13 = m->m12 * m->m23 * m->m44 -
          m->m12 * m->m43 * m->m24 -
          m->m13 * m->m22 * m->m44 +
          m->m13 * m->m42 * m->m24 +
          m->m14 * m->m22 * m->m43 -
          m->m14 * m->m42 * m->m23;

    inv.m14 = -m->m12 * m->m23 * m->m34 +
          m->m12 * m->m33 * m->m24 +
          m->m13 * m->m22 * m->m34 -
          m->m13 * m->m32 * m->m24 -
          m->m14 * m->m22 * m->m33 +
          m->m14 * m->m32 * m->m23;

    inv.m21 = -m->m21 * m->m33 * m->m44 +
          m->m21 * m->m43 * m->m34 +
          m->m23 * m->m31 * m->m44 -
          m->m23 * m->m41 * m->m34 -
          m->m24 * m->m31 * m->m43 +
          m->m24 * m->m41 * m->m33;

    inv.m22 = m->m11 * m->m33 * m->m44 -
          m->m11 * m->m43 * m->m34 -
          m->m13 * m->m31 * m->m44 +
          m->m13 * m->m41 * m->m34 +
          m->m14 * m->m31 * m->m43 -
          m->m14 * m->m41 * m->m33;

    inv.m23 = -m->m11 * m->m23 * m->m44 +
          m->m11 * m->m43 * m->m24 +
          m->m13 * m->m21 * m->m44 -
          m->m13 * m->m41 * m->m24 -
          m->m14 * m->m21 * m->m43 +
          m->m14 * m->m41 * m->m23;

    inv.m24 = m->m11 * m->m23 * m->m34 -
          m->m11 * m->m33 * m->m24 -
          m->m13 * m->m21 * m->m34 +
          m->m13 * m->m31 * m->m24 +
          m->m14 * m->m21 * m->m33 -
          m->m14 * m->m31 * m->m23;

    inv.m31 = m->m21 * m->m32 * m->m44 -
          m->m21 * m->m42 * m->m34 -
          m->m22 * m->m31 * m->m44 +
          m->m22 * m->m41 * m->m34 +
          m->m24 * m->m31 * m->m42 -
          m->m24 * m->m41 * m->m32;

    inv.m32 = -m->m11 * m->m32 * m->m44 +
          m->m11 * m->m42 * m->m34 +
          m->m12 * m->m31 * m->m44 -
          m->m12 * m->m41 * m->m34 -
          m->m14 * m->m31 * m->m42 +
          m->m14 * m->m41 * m->m32;

    inv.m33 = m->m11 * m->m22 * m->m44 -
          m->m11 * m->m42 * m->m24 -
          m->m12 * m->m21 * m->m44 +
          m->m12 * m->m41 * m->m24 +
          m->m14 * m->m21 * m->m42 -
          m->m14 * m->m41 * m->m22;

    inv.m34 = -m->m11 * m->m22 * m->m34 +
          m->m11 * m->m32 * m->m24 +
          m->m12 * m->m21 * m->m34 -
          m->m12 * m->m31 * m->m24 -
          m->m14 * m->m21 * m->m32 +
          m->m14 * m->m31 * m->m22;

    inv.m41 = -m->m21 * m->m32 * m->m43 +
          m->m21 * m->m42 * m->m33 +
          m->m22 * m->m31 * m->m43 -
          m->m22 * m->m41 * m->m33 -
          m->m23 * m->m31 * m->m42 +
          m->m23 * m->m41 * m->m32;

    inv.m42 = m->m11 * m->m32 * m->m43 -
          m->m11 * m->m42 * m->m33 -
          m->m12 * m->m31 * m->m43 +
          m->m12 * m->m41 * m->m33 +
          m->m13 * m->m31 * m->m42 -
          m->m13 * m->m41 * m->m32;

    inv.m43 = -m->m11 * m->m22 * m->m43 +
          m->m11 * m->m42 * m->m23 +
          m->m12 * m->m21 * m->m43 -
          m->m12 * m->m41 * m->m23 -
          m->m13 * m->m21 * m->m42 +
          m->m13 * m->m41 * m->m22;

    inv.m44 = m->m11 * m->m22 * m->m33 -
          m->m11 * m->m32 * m->m23 -
          m->m12 * m->m21 * m->m33 +
          m->m12 * m->m31 * m->m23 +
          m->m13 * m->m21 * m->m32 -
          m->m13 * m->m31 * m->m22;

    det = m->m11 * inv.m11 + m->m21 * inv.m12 + m->m31 * inv.m13 + m->m41 * inv.m14;

    /* Matrix Can Not Be Inverted (No Possible Inverse) if det == 0*/
    if (det != 0)
        det = 1.0 / det;

    result->m11 = inv.m11 * det;
    result->m21 = inv.m21 * det;
    result->m31 = inv.m31 * det;
    result->m41 = inv.m41 * det;

    result->m12 = inv.m12 * det;
    result->m22 = inv.m22 * det;
    result->m32 = inv.m32 * det;
    result->m42 = inv.m42 * det;

    result->m13 = inv.m13 * det;
    result->m23 = inv.m23 * det;
    result->m33 = inv.m33 * det;
    result->m43 = inv.m43 * det;

    result->m14 = inv.m14 * det;
    result->m24 = inv.m24 * det;
    result->m34 = inv.m34 * det;
    result->m44 = inv.m44 * det;
}

struct matrix matrix_inverse(struct matrix m)
{
    struct matrix result;
    pmatrix_inverse(&m, &result);
    return (result);
}
felselva commented 6 years ago

Thanks, will add soon!

Update: I added your name in CONTRIBUTORS.md.