estoy buscando una implementación de código de ejemplo sobre cómo invertir una matriz 4x4. Sé que hay eleminiación gaussiana, descomposición de LU, etc., pero en lugar de verlos en detalle, solo estoy buscando el código para hacerlo.invertir una matriz 4x4
idioma idealmente C++, los datos están disponibles en una matriz de 16 flotadores en orden cloumn-major.
gracias!
aquí:
bool gluInvertMatrix(const double m[16], double invOut[16])
{
double inv[16], det;
int i;
inv[0] = m[5] * m[10] * m[15] -
m[5] * m[11] * m[14] -
m[9] * m[6] * m[15] +
m[9] * m[7] * m[14] +
m[13] * m[6] * m[11] -
m[13] * m[7] * m[10];
inv[4] = -m[4] * m[10] * m[15] +
m[4] * m[11] * m[14] +
m[8] * m[6] * m[15] -
m[8] * m[7] * m[14] -
m[12] * m[6] * m[11] +
m[12] * m[7] * m[10];
inv[8] = m[4] * m[9] * m[15] -
m[4] * m[11] * m[13] -
m[8] * m[5] * m[15] +
m[8] * m[7] * m[13] +
m[12] * m[5] * m[11] -
m[12] * m[7] * m[9];
inv[12] = -m[4] * m[9] * m[14] +
m[4] * m[10] * m[13] +
m[8] * m[5] * m[14] -
m[8] * m[6] * m[13] -
m[12] * m[5] * m[10] +
m[12] * m[6] * m[9];
inv[1] = -m[1] * m[10] * m[15] +
m[1] * m[11] * m[14] +
m[9] * m[2] * m[15] -
m[9] * m[3] * m[14] -
m[13] * m[2] * m[11] +
m[13] * m[3] * m[10];
inv[5] = m[0] * m[10] * m[15] -
m[0] * m[11] * m[14] -
m[8] * m[2] * m[15] +
m[8] * m[3] * m[14] +
m[12] * m[2] * m[11] -
m[12] * m[3] * m[10];
inv[9] = -m[0] * m[9] * m[15] +
m[0] * m[11] * m[13] +
m[8] * m[1] * m[15] -
m[8] * m[3] * m[13] -
m[12] * m[1] * m[11] +
m[12] * m[3] * m[9];
inv[13] = m[0] * m[9] * m[14] -
m[0] * m[10] * m[13] -
m[8] * m[1] * m[14] +
m[8] * m[2] * m[13] +
m[12] * m[1] * m[10] -
m[12] * m[2] * m[9];
inv[2] = m[1] * m[6] * m[15] -
m[1] * m[7] * m[14] -
m[5] * m[2] * m[15] +
m[5] * m[3] * m[14] +
m[13] * m[2] * m[7] -
m[13] * m[3] * m[6];
inv[6] = -m[0] * m[6] * m[15] +
m[0] * m[7] * m[14] +
m[4] * m[2] * m[15] -
m[4] * m[3] * m[14] -
m[12] * m[2] * m[7] +
m[12] * m[3] * m[6];
inv[10] = m[0] * m[5] * m[15] -
m[0] * m[7] * m[13] -
m[4] * m[1] * m[15] +
m[4] * m[3] * m[13] +
m[12] * m[1] * m[7] -
m[12] * m[3] * m[5];
inv[14] = -m[0] * m[5] * m[14] +
m[0] * m[6] * m[13] +
m[4] * m[1] * m[14] -
m[4] * m[2] * m[13] -
m[12] * m[1] * m[6] +
m[12] * m[2] * m[5];
inv[3] = -m[1] * m[6] * m[11] +
m[1] * m[7] * m[10] +
m[5] * m[2] * m[11] -
m[5] * m[3] * m[10] -
m[9] * m[2] * m[7] +
m[9] * m[3] * m[6];
inv[7] = m[0] * m[6] * m[11] -
m[0] * m[7] * m[10] -
m[4] * m[2] * m[11] +
m[4] * m[3] * m[10] +
m[8] * m[2] * m[7] -
m[8] * m[3] * m[6];
inv[11] = -m[0] * m[5] * m[11] +
m[0] * m[7] * m[9] +
m[4] * m[1] * m[11] -
m[4] * m[3] * m[9] -
m[8] * m[1] * m[7] +
m[8] * m[3] * m[5];
inv[15] = m[0] * m[5] * m[10] -
m[0] * m[6] * m[9] -
m[4] * m[1] * m[10] +
m[4] * m[2] * m[9] +
m[8] * m[1] * m[6] -
m[8] * m[2] * m[5];
det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
if (det == 0)
return false;
det = 1.0/det;
for (i = 0; i < 16; i++)
invOut[i] = inv[i] * det;
return true;
}
Este fue levantada de MESA implementación de la biblioteca GLU.
Probablemente no lo quieras de otra manera. – shoosh
Sí, lo haría. Los compiladores son perfectamente capaces de desenrollar bucles, especialmente cuando se los indica. – Imagist
Incluso si no fuera así, podría ser el momento de invertir en la metaprogramación (aquí hay lugar para el ciclo). Resolví un compilador que no desenrollaba bucles triviales usando m4. –
DE MATRIZ 3x3
cambiar el código de acuerdo a sus necesidades
http://www.dreamincode.net/code/snippet1156.htm
Actualización:
Sí ... El ir de 3x3 a 4x4 parece una gran diferencia ... esta respuesta no es correcta para esto.
No hay necesidad de gritar. – jason
ir de 3x3 a 4x4 parece una gran diferencia – clamp
Para aclarar el punto de Matt, ese código menciona que calcula el determinante como parte del algoritmo para encontrar el inverso. Existen reglas simplificadas para calcular los determinantes de las matrices 2x2 y 3x3 que no se aplican a las matrices más grandes. – las3rjock
Puede usar el GNU Scientific Library o buscar el código en él.
Editar: Parece que quiere la sección Linear Algebra.
de hecho miré la matriz struct de gsl, pero no parece tener una función para determinante o inversión. – clamp
Aquí hay una biblioteca pequeña (solo un encabezado) C++ vector math (orientada a la programación 3D). Si lo usa, tenga en cuenta que el diseño de sus matrices en la memoria se invierte en comparación con lo que espera OpenGL, tuve tiempo divertido averiguarlo ...
Si necesita una biblioteca de C++ matriz con una gran cantidad de funciones, echar un vistazo a la biblioteca Eigen - http://eigen.tuxfamily.org
I 'enrollé' la aplicación MESA (también escribió un par de pruebas de unidad para asegurarse de que realmente funciona)
aquí:
float invf(int i,int j,const float* m){
int o = 2+(j-i);
i += 4+o;
j += 4-o;
#define e(a,b) m[ ((j+b)%4)*4 + ((i+a)%4) ]
float inv =
+ e(+1,-1)*e(+0,+0)*e(-1,+1)
+ e(+1,+1)*e(+0,-1)*e(-1,+0)
+ e(-1,-1)*e(+1,+0)*e(+0,+1)
- e(-1,-1)*e(+0,+0)*e(+1,+1)
- e(-1,+1)*e(+0,-1)*e(+1,+0)
- e(+1,-1)*e(-1,+0)*e(+0,+1);
return (o%2)?inv : -inv;
#undef e
}
bool inverseMatrix4x4(const float *m, float *out)
{
float inv[16];
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
inv[j*4+i] = invf(i,j,m);
double D = 0;
for(int k=0;k<4;k++) D += m[k] * inv[k*4];
if (D == 0) return false;
D = 1.0/D;
for (int i = 0; i < 16; i++)
out[i] = inv[i] * D;
return true;
}
escribí un poco sobre esto y mostrar el patrón de factores positivos/negativos on my blog.
Según lo sugerido por @LiraNuna, en muchas plataformas existen versiones aceleradas por hardware de tales rutinas, por lo que me complace tener una 'versión de copia de seguridad' que sea legible y concisa.
Nota: esto puede funcionar 3.5 veces más lento o peor que la implementación de MESA. Puede cambiar el patrón de factores para eliminar algunas adiciones, etc., pero perderá legibilidad y aún así no será muy rápido.
Puede hacerlo más rápido según este blog.
#define SUBP(i,j) input[i][j]
#define SUBQ(i,j) input[i][2+j]
#define SUBR(i,j) input[2+i][j]
#define SUBS(i,j) input[2+i][2+j]
#define OUTP(i,j) output[i][j]
#define OUTQ(i,j) output[i][2+j]
#define OUTR(i,j) output[2+i][j]
#define OUTS(i,j) output[2+i][2+j]
#define INVP(i,j) invP[i][j]
#define INVPQ(i,j) invPQ[i][j]
#define RINVP(i,j) RinvP[i][j]
#define INVPQ(i,j) invPQ[i][j]
#define RINVPQ(i,j) RinvPQ[i][j]
#define INVPQR(i,j) invPQR[i][j]
#define INVS(i,j) invS[i][j]
#define MULTI(MAT1, MAT2, MAT3) \
MAT3(0,0)=MAT1(0,0)*MAT2(0,0) + MAT1(0,1)*MAT2(1,0); \
MAT3(0,1)=MAT1(0,0)*MAT2(0,1) + MAT1(0,1)*MAT2(1,1); \
MAT3(1,0)=MAT1(1,0)*MAT2(0,0) + MAT1(1,1)*MAT2(1,0); \
MAT3(1,1)=MAT1(1,0)*MAT2(0,1) + MAT1(1,1)*MAT2(1,1);
#define INV(MAT1, MAT2) \
_det = 1.0/(MAT1(0,0) * MAT1(1,1) - MAT1(0,1) * MAT1(1,0)); \
MAT2(0,0) = MAT1(1,1) * _det; \
MAT2(1,1) = MAT1(0,0) * _det; \
MAT2(0,1) = -MAT1(0,1) * _det; \
MAT2(1,0) = -MAT1(1,0) * _det; \
#define SUBTRACT(MAT1, MAT2, MAT3) \
MAT3(0,0)=MAT1(0,0) - MAT2(0,0); \
MAT3(0,1)=MAT1(0,1) - MAT2(0,1); \
MAT3(1,0)=MAT1(1,0) - MAT2(1,0); \
MAT3(1,1)=MAT1(1,1) - MAT2(1,1);
#define NEGATIVE(MAT) \
MAT(0,0)=-MAT(0,0); \
MAT(0,1)=-MAT(0,1); \
MAT(1,0)=-MAT(1,0); \
MAT(1,1)=-MAT(1,1);
void getInvertMatrix(complex<double> input[4][4], complex<double> output[4][4]) {
complex<double> _det;
complex<double> invP[2][2];
complex<double> invPQ[2][2];
complex<double> RinvP[2][2];
complex<double> RinvPQ[2][2];
complex<double> invPQR[2][2];
complex<double> invS[2][2];
INV(SUBP, INVP);
MULTI(SUBR, INVP, RINVP);
MULTI(INVP, SUBQ, INVPQ);
MULTI(RINVP, SUBQ, RINVPQ);
SUBTRACT(SUBS, RINVPQ, INVS);
INV(INVS, OUTS);
NEGATIVE(OUTS);
MULTI(OUTS, RINVP, OUTR);
MULTI(INVPQ, OUTS, OUTQ);
MULTI(INVPQ, OUTR, INVPQR);
SUBTRACT(INVP, INVPQR, OUTP);
}
Esto no es una implementación completa porque P puede que no sea invertible, pero se puede combinar este código con la aplicación de MESA para obtener un mejor rendimiento.
Si cualquiera que busque un código más costumized y "fácil de leer", entonces me consiguió este
var A2323 = m.m22 * m.m33 - m.m23 * m.m32 ;
var A1323 = m.m21 * m.m33 - m.m23 * m.m31 ;
var A1223 = m.m21 * m.m32 - m.m22 * m.m31 ;
var A0323 = m.m20 * m.m33 - m.m23 * m.m30 ;
var A0223 = m.m20 * m.m32 - m.m22 * m.m30 ;
var A= m.m20 * m.m31 - m.m21 * m.m30 ;
var A2313 = m.m12 * m.m33 - m.m13 * m.m32 ;
var A1313 = m.m11 * m.m33 - m.m13 * m.m31 ;
var A1213 = m.m11 * m.m32 - m.m12 * m.m31 ;
var A2312 = m.m12 * m.m23 - m.m13 * m.m22 ;
var A1312 = m.m11 * m.m23 - m.m13 * m.m21 ;
var A1212 = m.m11 * m.m22 - m.m12 * m.m21 ;
var A0313 = m.m10 * m.m33 - m.m13 * m.m30 ;
var A0213 = m.m10 * m.m32 - m.m12 * m.m30 ;
var A0312 = m.m10 * m.m23 - m.m13 * m.m20 ;
var A0212 = m.m10 * m.m22 - m.m12 * m.m20 ;
var A0113 = m.m10 * m.m31 - m.m11 * m.m30 ;
var A0112 = m.m10 * m.m21 - m.m11 * m.m20 ;
var det = m.m00 * (m.m11 * A2323 - m.m12 * A1323 + m.m13 * A1223)
- m.m01 * (m.m10 * A2323 - m.m12 * A0323 + m.m13 * A0223)
+ m.m02 * (m.m10 * A1323 - m.m11 * A0323 + m.m13 * A0123)
- m.m03 * (m.m10 * A1223 - m.m11 * A0223 + m.m12 * A0123) ;
det = 1/det;
return new Matrix4x4() {
m00 = det * (m.m11 * A2323 - m.m12 * A1323 + m.m13 * A1223),
m01 = det * - (m.m01 * A2323 - m.m02 * A1323 + m.m03 * A1223),
m02 = det * (m.m01 * A2313 - m.m02 * A1313 + m.m03 * A1213),
m03 = det * - (m.m01 * A2312 - m.m02 * A1312 + m.m03 * A1212),
m10 = det * - (m.m10 * A2323 - m.m12 * A0323 + m.m13 * A0223),
m11 = det * (m.m00 * A2323 - m.m02 * A0323 + m.m03 * A0223),
m12 = det * - (m.m00 * A2313 - m.m02 * A0313 + m.m03 * A0213),
m13 = det * (m.m00 * A2312 - m.m02 * A0312 + m.m03 * A0212),
m20 = det * (m.m10 * A1323 - m.m11 * A0323 + m.m13 * A0123),
m21 = det * - (m.m00 * A1323 - m.m01 * A0323 + m.m03 * A0123),
m22 = det * (m.m00 * A1313 - m.m01 * A0313 + m.m03 * A0113),
m23 = det * - (m.m00 * A1312 - m.m01 * A0312 + m.m03 * A0112),
m30 = det * - (m.m10 * A1223 - m.m11 * A0223 + m.m12 * A0123),
m31 = det * (m.m00 * A1223 - m.m01 * A0223 + m.m02 * A0123),
m32 = det * - (m.m00 * A1213 - m.m01 * A0213 + m.m02 * A0113),
m33 = det * (m.m00 * A1212 - m.m01 * A0212 + m.m02 * A0112),
};
no escribo el código, pero mi programa lo hice. Hice un pequeño programa para hacer un programa que calcule el determinante e inverso de cualquier matriz N.
Lo hago porque una vez en el pasado necesito un código que invierte la matriz de 5x5, pero nadie en la tierra ha hecho esto, así que hice uno.
Eche un vistazo sobre el programa here.
¿Es esta tarea? Si no es así (por ejemplo, solo intentas resolver Ax = b), intentar calcular explícitamente un inverso puede no ser lo que quieres hacer. –
no es tarea. es para un proyecto personal. y no quiero "perder" el tiempo aprendiendo la inversión de matriz para 4x4, lo que parece bastante complicado en comparación con 3x3 – clamp
No creo que esta sea una pregunta estúpida que merece -1 puntaje. – stribika