\begin{equation}\label{eq:1}A^{-1} = \frac {1}{det(A)}\cdot C_{A}^T\end{equation}
dove
\begin{equation}\label{eq:2}det(A)=\sum_{j=0}^{n-1}[(-1)^{i+j}\cdot a_{i,j}\cdot det(A_{i,j})]\end{equation}
con $0 \leq i < n$ fissato e costante, con $det(A_{i,j})$ minore di $A$ (ottenuto eliminando la $i$-esima riga e la $j$-esima colonna) e con $a_{i,j}$ elemento della matrice $A$ sulla $i$-esima riga e $j$-esima colonna.
$C_{A}^T$ in \eqref{eq:1}, invece, è la trasposta della matrice dei cofattori di $A$, chiamata anche matrice aggiunta. Ogni elemento di $C_{A}^T$ si può definire quindi con
\begin{equation}\label{eq:3}c_{i,j}^{T} = (-1)^{i+j} \cdot det(A_{j,i})\end{equation}
dove $det(A_{j,i})$ è il minore di $A$ ottenuto eliminando la $j$-esima riga e la $i$-esima colonna (si noti l'inversione degli indici per selezionare la trasposta; senza tale inversione si otterrebbe $C_A$, la matrice dei cofattori di $A$).
Guardando le prime tre equazioni si può definire e calcolare ogni elemento dell'inversa di $A$ come
\begin{equation}\label{eq:4}a_{i,j}^{-1} = \frac {c_{i,j}^{T}}{det(A)} \end{equation}
dove $a_{i,j}^{-1}$ è un elemento di $A^{-1}$ e $c_{i,j}^{T}$ è un elemento di $C_{A}^T$.
In Computer Graphics si ha a che fare perlopiù come matrici $4\times4$ ed in questo caso un primo abbozzo di codice molto specializzato potrebbe essere
/*Per facilitare la stesura, la lettura e la comprensione del codice si è volutamente deciso di utilizzare un mix di C/C++ cercando di restare il più possibile all'interno del sottoinsieme C del C++ e di utilizzare solo quei costrutti e funzionalità di base del C++ che permettono una maggiore semplificazione.*/ #include <stdio.h> #include <math.h> void inverse(float*, float); void adjoint(float*); float det4x4(float*); float det3x3(float, float, float, float, float, float, float, float, float); float det2x2(float, float, float, float); void print_matrix(float*, int); int main(int argc, char *argv[], char *envp[]) { float m4x4[4][4] = { {3, -0.1, -0.2, 0}, {0.1, 7, -0.3, 0}, {0.3, -0.2, 10, 0}, {0, 0, 0, 1} }; float adj4x4[4][4]; // Se si vuole conservare la matrice m4x4 è necessario copiare elementi in matrice di output adj4x4 for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) adj4x4[i][j] = m4x4[i][j]; } inverse(&adj4x4[0][0], 1e-6); print_matrix(&adj4x4[0][0], 4); return 0; } void inverse(float *a, float tol) { float det; det = det4x4(a); adjoint(a); if (fabs(det) < tol) fprintf(stderr, "Errore di tolleranza; prababile matrice singolare"); for (int i = 0; i < 4; ++i) for (int j = 0; j < 4; ++j) a[4 * i + j] = a[4 * i + j] / det; // vedi (4) } void adjoint(float* a) { // Da un nome ad ogni elemento della matrice per convenienza float a0, a1, a2, a3; float b0, b1, b2, b3; float c0, c1, c2, c3; float d0, d1, d2, d3; a0 = a[4 * 0 + 0]; a1 = a[4 * 0 + 1]; a2 = a[4 * 0 + 2]; a3 = a[4 * 0 + 3]; b0 = a[4 * 1 + 0]; b1 = a[4 * 1 + 1]; b2 = a[4 * 1 + 2]; b3 = a[4 * 1 + 3]; c0 = a[4 * 2 + 0]; c1 = a[4 * 2 + 1]; c2 = a[4 * 2 + 2]; c3 = a[4 * 2 + 3]; d0 = a[4 * 3 + 0]; d1 = a[4 * 3 + 1]; d2 = a[4 * 3 + 2]; d3 = a[4 * 3 + 3]; // Calcola elementi di matrice aggiunta; vedi eq. (3) // Attenzione: sovrascrive matrice di input // Scambia righe e colonne per ottenere matrice aggiunta a[4 * 0 + 0] = det3x3(b1, b2, b3, c1, c2, c3, d1, d2, d3); a[4 * 0 + 1] = -det3x3(a1, a2, a3, c1, c2, c3, d1, d2, d3); a[4 * 0 + 2] = det3x3(a1, a2, a3, b1, b2, b3, d1, d2, d3); a[4 * 0 + 3] = -det3x3(a1, a2, a3, b1, b2, b3, c1, c2, c3); a[4 * 1 + 0] = -det3x3(b0, b2, b3, c0, c2, c3, d0, d2, d3); a[4 * 1 + 1] = det3x3(a0, a2, a3, c0, c2, c3, d0, d2, d3); a[4 * 1 + 2] = -det3x3(a0, a2, a3, b0, b2, b3, d0, d2, d3); a[4 * 1 + 3] = det3x3(a0, a2, a3, b0, b2, b3, c0, c2, c3); a[4 * 2 + 0] = det3x3(b0, b1, b3, c0, c1, c3, d0, d1, d3); a[4 * 2 + 1] = -det3x3(a0, a1, a3, c0, c1, c3, d0, d1, d3); a[4 * 2 + 2] = det3x3(a0, a1, a3, b0, b1, b3, d0, d1, d3); a[4 * 2 + 3] = -det3x3(a0, a1, a3, b0, b1, b3, c0, c1, c3); a[4 * 3 + 0] = -det3x3(b0, b1, b2, c0, c1, c2, d0, d1, d2); a[4 * 3 + 1] = det3x3(a0, a1, a2, c0, c1, c2, d0, d1, d2); a[4 * 3 + 2] = -det3x3(a0, a1, a2, b0, b1, b2, d0, d1, d2); a[4 * 3 + 3] = det3x3(a0, a1, a2, b0, b1, b2, c0, c1, c2); } float det4x4(float *a) { float det; // Da un nome ad ogni elemento della matrice per convenienza float a0, a1, a2, a3; float b0, b1, b2, b3; float c0, c1, c2, c3; float d0, d1, d2, d3; a0 = a[4 * 0 + 0]; a1 = a[4 * 0 + 1]; a2 = a[4 * 0 + 2]; a3 = a[4 * 0 + 3]; b0 = a[4 * 1 + 0]; b1 = a[4 * 1 + 1]; b2 = a[4 * 1 + 2]; b3 = a[4 * 1 + 3]; c0 = a[4 * 2 + 0]; c1 = a[4 * 2 + 1]; c2 = a[4 * 2 + 2]; c3 = a[4 * 2 + 3]; d0 = a[4 * 3 + 0]; d1 = a[4 * 3 + 1]; d2 = a[4 * 3 + 2]; d3 = a[4 * 3 + 3]; // vedi eq. (2) det = a0 * det3x3(b1, b2, b3, c1, c2, c3, d1, d2, d3) - a1 * det3x3(b0, b2, b3, c0, c2, c3, d0, d2, d3) + a2 * det3x3(b0, b1, b3, c0, c1, c3, d0, d1, d3) - a3 * det3x3(b0, b1, b2, c0, c1, c2, d0, d1, d2); return det; } float det3x3(float a0, float a1, float a2, float b0, float b1, float b2, float c0, float c1, float c2) { float det; // vedi eq. (2) det = a0 * det2x2(b1, b2, c1, c2) - a1 * det2x2(b0, b2, c0, c2) + a2 * det2x2(b0, b1, c0, c1); return det; } float det2x2(float b0, float b1, float c0, float c1) { // vedi eq. (2) return b0 * c1 - c0 * b1; } void print_matrix(float *m, int n) { for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { printf("%f ", m[n*i + j]); } printf("\n"); } printf("\n"); }
Non il massimo dell'ottimizzazione ma funziona per qualsiasi matrice $4\times4$ ed è relativamente efficiente (volendo si possono eliminare le chiamate a funzione e mettere tutto in inverse(), oppure si possono rendere inline; per queste ed altre ottimizzazioni si vedano [3] e [4]).
Un reale miglioramento delle performance si può ottenere se si considerano solo matrici collegate a trasformazioni affini, che sono poi le più utilizzate in computer graphics. In questo caso l'ultima riga (in OpenGL, colonna in DirectX) sarà sempre $[0, 0, 0, 1]$ e l'intera matrice si può dividere a blocchi in questo modo
\begin{equation}\label{eq:5}\begin{bmatrix} \mathbf{A} & \mathbf{t}\\ \mathbf{0}^T & 1\\ \end{bmatrix} \quad\quad \mbox{in OpenGL} \end{equation}
\begin{equation}\label{eq:6}\begin{bmatrix} \mathbf{A} & \mathbf{0}\\ \mathbf{t}^T & 1\\ \end{bmatrix} \quad\quad \mbox{in DirectX} \end{equation}
dove $\mathbf{A}$ è la matrice superiore destra $3\times3$ che identifica rotazioni e scaling, $\mathbf{t}$ è il vettore colonna che identifica le traslazioni, $\mathbf{0}$ è il vettore colonna nullo e $1$ è una costante scalare.
Per quanto visto in [6], l'inversa di una matrice affine è definita come
\begin{equation}\label{eq:7} \begin{bmatrix} \mathbf{A} & \mathbf{t} \\ \mathbf{0}^T & 1 \\ \end{bmatrix}^{-1} = \begin{bmatrix} \mathbf{A^{-1}} & -\mathbf{A^{-1}}\mathbf{t} \\ \mathbf{0}^T & 1 \\ \end{bmatrix} \quad\quad \mbox{In OpenGL}\end{equation}
\[ \begin{bmatrix} \mathbf{A} & \mathbf{0} \\ \mathbf{t}^T & 1 \\ \end{bmatrix}^{-1} = \begin{bmatrix} \mathbf{A^{-1}} & \mathbf{0} \\ (-\mathbf{A^{-1}}\mathbf{t})^T & 1 \\ \end{bmatrix} \quad\quad \mbox{In DirectX} \]
In pratica è sufficiente invertire la matrice $3\times3$ superiore sinistra per trovare la relativa matrice superiore sinistra dell'inversa. In seguito basterà moltiplicarla per $\mathbf{t}$ per trovare i primi 3 elementi della quarta colonna (in OpenGL, riga in DirectX). Infine non resta che impostare l'ultima riga (in OpenGL, colonna in DirectX) a $[0, 0, 0, 1]$ poichè questa non cambia nelle trasformazioni affini.
void inverse(float *a, float tol) { float det; // Da un nome ad ogni elemento della matrice per convenienza float a0, a1, a2, a3; float b0, b1, b2, b3; float c0, c1, c2, c3; float d0, d1, d2, d3; a0 = a[4 * 0 + 0]; a1 = a[4 * 0 + 1]; a2 = a[4 * 0 + 2]; a3 = a[4 * 0 + 3]; b0 = a[4 * 1 + 0]; b1 = a[4 * 1 + 1]; b2 = a[4 * 1 + 2]; b3 = a[4 * 1 + 3]; c0 = a[4 * 2 + 0]; c1 = a[4 * 2 + 1]; c2 = a[4 * 2 + 2]; c3 = a[4 * 2 + 3]; d0 = a[4 * 3 + 0]; d1 = a[4 * 3 + 1]; d2 = a[4 * 3 + 2]; d3 = a[4 * 3 + 3]; // Determinante matrice superiore sinistra det = det3x3(a0, a1, a2, b0, b1, b2, c0, c1, c2); if (fabs(det) < tol) fprintf(stderr, "Errore di tolleranza; prababile matrice singolare"); // Inversa di matrice superiore sinistra; A^1 in eq. (7) // Attenzione: sovrascrive matrice di input a[4 * 0 + 0] = det2x2(b1, b2, c1, c2) / det; a[4 * 0 + 1] = -(det2x2(a1, a2, c1, c2) / det); a[4 * 0 + 2] = det2x2(a1, a2, b1, b2) / det; a[4 * 1 + 0] = -(det2x2(b0, b2, c0, c2) / det); a[4 * 1 + 1] = det2x2(a0, a2, c0, c2) / det; a[4 * 1 + 2] = -(det2x2(a0, a2, b0, b2) / det); a[4 * 2 + 0] = det2x2(b0, b1, c0, c1) / det; a[4 * 2 + 1] = -(det2x2(a0, a1, c0, c1) / det); a[4 * 2 + 2] = det2x2(a0, a1, b0, b1) / det; // -A^1*c in eq. (7) // Attenzione: sovrascrive matrice di input a[4 * 0 + 3] = -(a[4 * 0 + 0] * a[4 * 0 + 3] + a[4 * 0 + 1] * a[4 * 1 + 3] + a[4 * 0 + 2] * a[4 * 2 + 3]); a[4 * 1 + 3] = -(a[4 * 1 + 0] * a[4 * 0 + 3] + a[4 * 1 + 1] * a[4 * 1 + 3] + a[4 * 1 + 2] * a[4 * 2 + 3]); a[4 * 2 + 3] = -(a[4 * 2 + 0] * a[4 * 0 + 3] + a[4 * 2 + 1] * a[4 * 1 + 3] + a[4 * 2 + 2] * a[4 * 2 + 3]); // Ultima riga rimane 0, 0, 0, 1 // Attenzione: sovrascrive matrice di input a[4 * 3 + 0] = a[4 * 3 + 1] = a[4 * 3 + 2] = 0; a[4 * 3 + 3] = 1; } float det3x3(float a0, float a1, float a2, float b0, float b1, float b2, float c0, float c1, float c2) { float det; // vedi eq. (2) det = a0 * det2x2(b1, b2, c1, c2) - a1 * det2x2(b0, b2, c0, c2) + a2 * det2x2(b0, b1, c0, c1); return det; } float det2x2(float b0, float b1, float c0, float c1) { // vedi eq. (2) return b0 * c1 - c0 * b1; }
Infine, questo è il confronto tra il risultato del nostro codice e quello restituito da Wolfram|Alpha
Riferimenti
[1] Graphics Gems I - Glassner
[2] Graphics Gems II - Arvo
[3] Foundations of Game Engine Development Vol.1 - Lengyel
[4] Game Physics Engine Development - Millington
[5] https://it.wikipedia.org/wiki/Matrice_invertibile
[6] Trasformazioni Lineari ed Affini in Computer Graphics
Nessun commento:
Posta un commento