3D-Programmierung (Mathematik) Teil 07: inverse Matrix versus transponierte Matrix

Teil 7 dieser Artikelreihe ist für alle diejenigen Programmierer gedacht, die bereits an einer eigenen Mathematik-Bibliothek arbeiten bzw. sich demnächst damit befassen wollen. Im Folgenden werden beispielhafte Implementierungen zweier wichtiger Matrizenfunktionen vorgestellt. Die erste Funktion dient zur Berechnung einer transponierten Matrix, die zweite zur Berechnung einer inversen Matrix.


Hinweis:

Mithilfe von transponierten Matrizen lassen sich Rotationstransformationen wieder rückgängig machen. Bewirkt eine Rotationsmatrix beispielsweise eine Drehung von 10° um die x-Achse, so bewirkt die transponierte Rotationsmatrix eine Drehung von –10° um dieselbe Achse. Eine invertierte Rotationsmatrix hätte zwar denselben Effekt, jedoch ist der Berechnungsaufwand für eine inverse Matrix deutlich höher. Sollen jedoch Rotations-Translations-Transformationen rückgängig gemacht werden, dann kommt man um die Verwendung der inversen Matrix nicht herum.

transponierte Matrix

INLINE D3DXMATRIX* D3DXMatrixTranspose(D3DXMATRIX *pOut,
                                       CONST D3DXMATRIX *pM)
{

    int i,j;

    for(i=0; i<4; i++)
    {
        for(j=0; j<4; j++)
        {
            pOut->m[i][j] = pM->m[j][i];
        }
    }

    return pOut;
}

inverse Matrix

INLINE D3DXMATRIX* D3DXMatrixInverse(D3DXMATRIX *pOut,
                                     FLOAT *pDeterminant,
                                     CONST D3DXMATRIX *pM)

{
    int a, i, j;
    D3DXMATRIX out;
    D3DXVECTOR4 v, vec[3];
    FLOAT det;

    det = D3DXMatrixDeterminant(pM);

    if(!det )
        return NULL;

    if(pDeterminant)
        *pDeterminant = det;

    for(i=0; i<4; i++)
    {
        for(j=0; j<4; j++)
        {
            if(j != i)
            {
                a = j;

                if(j > i)
                    a = a-1;

                vec[a].x = pM->m[j][0];
                vec[a].y = pM->m[j][1];
                vec[a].z = pM->m[j][2];
                vec[a].w = pM->m[j][3];
            }
        }

        D3DXVec4Cross(&v, &vec[0], &vec[1], &vec[2]);

        // pow(-1.0f, i) = (-1.0f) hoch i
        out.m[0][i] = pow(-1.0f, i) * v.x / det;
        out.m[1][i] = pow(-1.0f, i) * v.y / det;
        out.m[2][i] = pow(-1.0f, i) * v.z / det;
        out.m[3][i] = pow(-1.0f, i) * v.w / det;
    }

    *pOut = out;
    return pOut;
}

Weitere wichtige Funktionen zur Berechnung einer inversen Matrix

Determinate einer Matrix

INLINE FLOAT D3DXMatrixDeterminant(CONST D3DXMATRIX *pM)
{
    D3DXVECTOR4 minor, v1, v2, v3;
    FLOAT det;

    v1.x = pM->m[0][0];
    v1.y = pM->m[1][0];
    v1.z = pM->m[2][0];
    v1.w = pM->m[3][0];

    v2.x = pM->m[0][1];
    v2.y = pM->m[1][1];
    v2.z = pM->m[2][1];
    v2.w = pM->m[3][1];

    v3.x = pM->m[0][2];
    v3.y = pM->m[1][2];
    v3.z = pM->m[2][2];
    v3.w = pM->m[3][2];

    D3DXVec4Cross(&minor,&v1,&v2,&v3);

    det = -(pM->m[0][3] * minor.x +
            pM->m[1][3] * minor.y +
            pM->m[2][3] * minor.z +
            pM->m[3][3] * minor.w);

    return det;
}

Kreuzprodukt dreier vierdimensionaler Vektoren

INLINE D3DXVECTOR4* D3DXVec4Cross(D3DXVECTOR4 *pOut,
                                  CONST D3DXVECTOR4 *pV1,
                                  CONST D3DXVECTOR4 *pV2,
                                  CONST D3DXVECTOR4 *pV3)
{
    pOut->x = pV1->y * (pV2->z * pV3->w - pV3->z * pV2->w) -
              pV1->z * (pV2->y * pV3->w - pV3->y * pV2->w) +
              pV1->w * (pV2->y * pV3->z - pV2->z *pV3->y);

    pOut->y = -(pV1->x * (pV2->z * pV3->w - pV3->z * pV2->w) -
                pV1->z * (pV2->x * pV3->w - pV3->x * pV2->w) +
                pV1->w * (pV2->x * pV3->z - pV3->x * pV2->z));

    pOut->z = pV1->x * (pV2->y * pV3->w - pV3->y * pV2->w) -
              pV1->y * (pV2->x *pV3->w - pV3->x * pV2->w) +
              pV1->w * (pV2->x * pV3->y - pV3->x * pV2->y);

    pOut->w = -(pV1->x * (pV2->y * pV3->z - pV3->y * pV2->z) -
                pV1->y * (pV2->x * pV3->z - pV3->x * pV2->z) +
                pV1->z * (pV2->x * pV3->y - pV3->x * pV2->y));

    return pOut;
}