Matriz inversa e intersección de dos planos

[Findeton]

                                Hola,estoy intentando crear una serie de funciones que me sirvan para detectar la colisión entre polígonos, esferas y líneas y me he encontrado un problema:

Estoy intentando encontrar la intersección entre 2 planos, los cuales han sido calculados a partir de 2 polígonos, dándome así la ecuación de cada plano Ax +By + Cz +D = 0, siendo el vector normal de cada plano (x=A, y=B, z=C). Entonces, utilizando el método de Microsoft, que viene descrito en la web http://research.microsoft.com/~jckrumm/Int..._two_planes.htm , sólo necesito resolver un sistema de ecuaciones en forma de matrices.


Para ello, me fui a la web http://www.sc.ehu.es/sbweb/fisica/cursoJav...iz1/matriz1.htm donde explican que para saber la solución de A*S=B (donde A es una matriz conocida, B un vector conocido y S el vector solución), hay que calcular la inversa de A y multiplicarla por B, y te dará S. Entonces, como en esa misma web me viene el código para calcular la inversa... pensé que todo estaba solucionado... pasé el código a C++ y... pues sólo calcula bien la inversa de matrices 2x2, cuando el método descrito en la web es para nxn, donde n es cualquier número entero.

Así que posteo el código para que comprobéis si hay algún error, si podéis  :(.

Ésta es la función que calcula la matriz inversa:

Código [Seleccionar] Expand



Matriz matriz_inversa(Matriz d){

       int n=d.n;     //dimensión de la matriz

 Matriz a(n);

 a=d;

       Matriz b(n);    //matriz de los términos independientes

       Matriz c(n);    //matriz de las incógnitas



//matriz unidad

       for(int i=0; i<n; i++){

           b.x[i][i]=1.0;

       }

//transformación de la matriz y de los términos independientes

       for(int k=0; k<n-1; k++){

           for(int i=k+1; i<n; i++){

//términos independientes

               for(int s=0; s<n; s++){

                   b.x[i][s]-=a.x[i][k]*b.x[k][s]/a.x[k][k];

               }

//elementos de la matriz

               for(int j=k+1; j<n; j++){

                   a.x[i][j]-=a.x[i][k]*a.x[k][j]/a.x[k][k];

               }

           }

       }

//cálculo de las incógnitas, elementos de la matriz inversa

       for(int s=0; s<n; s++){

           c.x[n-1][s]=b.x[n-1][s]/a.x[n-1][n-1];

           for(int i=n-2; i>=0; i--){

               c.x[i][s]=b.x[i][s]/a.x[i][i];

               for(int k=n-1; k>i; k--){

                   c.x[i][s]-=a.x[i][k]*c.x[k][s]/a.x[i][i];

               }

           }

       }

       return c;

   }



Y aquí está la definicíon de la clase Matriz:

Código [Seleccionar] Expand



class Matriz {



public:

int n;         //dimensión

   double** x;     //el valor

   Matriz(int m); //constructor

  ~Matriz();     //destructor

  void mostrar();

   };





//Aquí definimos el constructor



Matriz::Matriz(int m) {

       n=m;

       x= new double*[n];

       for(int i=0; i<n; i++){

  x[i]= new double[n];

 }  

       for(i=0; i<n; i++){

           for(int j=0; j<n; j++){

               x[i][j]=0.0;

           }

       }

   }



//Aquí definimos el destructor



Matriz::~Matriz() {

   

                     

   }



void Matriz::identidad(){

       for(int i=0; i<n; i++){

           for(int j=0; j<n; j++){

               x[i][j]=0.0;

           }

       }

 for(i=0; i<n; i++)

  x[i][i]=1;

}

void Matriz::identidad(double m){

       for(int i=0; i<n; i++){

           for(int j=0; j<n; j++){

               x[i][j]=0.0;

           }

       }

 for(i=0; i<n; i++)

  x[i][i]=m;

}



void Matriz::mostrar(){



for(int i=0; i<n;i++){ //hay que poner '< n' y no '<= n' pq si no, se pasará, ya que las matrices empiezan en 0

 for(int k=0; k<n;k++){//hay que poner '< n' y no '<= n' pq si no, se pasará, ya que las matrices empiezan en 0

        cout << x[i][k] << " ";

 }

 cout << endl;

}

}






Si necesitáis más datos, por favor pedirlos, pero ya os digo que un código de ejemplo sería:

Código [Seleccionar] Expand



#include "iostream.h"



/*Poned aquí la definición de la matriz y sus funciones*/



int main(int argc, char* argv[])

{

Matriz b(3);



b.x[0][0]=4;

b.x[0][1]=1;

b.x[0][2]=1;

b.x[1][0]=2;

b.x[1][1]=3;

b.x[1][2]=1;

b.x[2][0]=2;

b.x[2][1]=3;

b.x[2][2]=1;



   

b=matriz_inversa(b);

b.mostrar();



return 0;

}                                

[Juan Mellado]

Hola Findeton.

El código parece que está correctamente portado. No lo he mirado a fondo, es que al probar me he dado cuenta que la matriz que estás utilizando no es invertible. Es decir, su determinante es cero.

De forma general: una matriz es invertible si y sólo si su determinante es distinto de cero.

Espero haberte ayudado.
Saludos a todos.

[Findeton]

                                Pues sí que me has ayudados, ahora lo que me gustaría saber es si una matriz en principio no invertible puede ser transformada en otra equivalente sí invertible. Pero ya buscaré en la web.                                

[Ryuchan]

                                Si el sistema de referencia que usas es ortonormal (que seguro que lo es (ejes perpendiculares 2 a 2 y de modulo 1)), la inversa de tu matriz es la transpuesta.
Asi de simple :)                                

[Findeton]

                                Entonces en toda matriz basada en un ssitema de referencia ortogonal se cumple que la transpuesta de la misma es igual a su inversa?                                

[Mars Attacks]

                                Inversa de matriz con bases ortonormales igual a matriz traspuesta. Vamos, que sí.                                

[Findeton]

                                Pues sí, finalmente conseguí crearme mi propia función para detectar la colisión entre 2 planos. La pongo a continuación por si a alguien le ayuda:

Código [Seleccionar] Expand





// Aquí declaramos las clases:



class Matriz {



public:

int n;         //dimensión

   double** x;     //el valor

   Matriz(int m); //constructor

  ~Matriz();     //destructor

  void mostrar();

   };



//Aquí definimos el constructor



Matriz::Matriz(int m) {

       n=m;

       x= new double*[n];

       for(int i=0; i<n; i++){

  x[i]= new double[n];

 }  

       for(i=0; i<n; i++){

           for(int j=0; j<n; j++){

               x[i][j]=0.0;

           }

       }

   }



//Aquí definimos el destructor



Matriz::~Matriz() {

   

                     

   }



void Matriz::mostrar(){



for(int i=0; i<n;i++){ //hay que poner '< n' y no '<= n' pq si no, se pasará, ya que las matrices empiezan en 0

 for(int k=0; k<n;k++){//hay que poner '< n' y no '<= n' pq si no, se pasará, ya que las matrices empiezan en 0

        cout << x[i][k] << " ";

 }

 cout << endl;

}

}



class Vector {

public:

 int n;      //dimensión

   double* x;

   Vector(int m);



~Vector();

void mostrar();



   

//otras funciones miembro

};





Vector::Vector(int m) {

       n=m;

       x=new double[n];

       for(int i=0; i<n; i++)//hay que poner '< n' y no '<= n' pq si no, se pasará, ya que las matrices empiezan en 0

           x[i]=0;

       

   }







Vector::~Vector() {

       

             

   }



void Vector::mostrar(){



for(int i=0; i<n;i++){

cout << x[i] << ", ";

}

}



class punto3d {

public:

double x, y, z;

punto3d();

};



punto3d::punto3d(){

x=y=z=0;

}



class vector3d {

public:

double x, y, z;

normalizar();

vector3d();

};



vector3d::vector3d(){

x=y=z=0;

}



vector3d::normalizar() {

double len = (double)sqrt(x*x + y*y + z*z);

x /= len;

y /= len;

z /= len;

}





class vlinea3d{

public:

punto3d p;

vector3d v;

};





class plano{

public:

 double A, B, C, D;

 void coeficientes(punto3d a, punto3d b, punto3d  c);

 vector3d v;

 void vector();

 plano();



};



plano::plano(){

A=B=C=D=0;

}



void plano::vector(){

v.x=A;

v.y=B;

v.z=C;



}

/// Coeficietes calcula los coeficientes del plano                                  ///

void plano::coeficientes(punto3d a, punto3d b, punto3d  c){

double A=0;

A=a.y*(b.z-c.z) + b.y*(c.z-a.z) + c.y*(a.z-b.z);

B=a.z*(b.x-c.x) + b.z*(c.x-a.x) + c.z*(a.x-b.x);

C=a.x*(b.y-c.y) + b.x*(c.y-a.y) + c.x*(a.y-b.y);

D=(a.x*(b.y*c.z-c.y*b.z)+b.x*(c.y*a.z-a.y*c.z)+c.x*(a.y*b.z-b.y*a.z));

}



class poligono3d{

public:

punto3d A,B,C;

plano plano;

               poligono3d();

};



poligono3d::poligono3d(){

A.x=A.y=A.z=B.x=B.y=B.z=C.x=C.y=C.z=0;

}





inline vector3d CrossProduct (vector3d A, vector3d B) {

vector3d C;

C.x = A.y*B.z - A.z*B.y;

C.y = A.z*B.x - A.x*B.z;

C.z = A.x*B.y - A.y*B.x;



return C;

}



inline Vector vector_producto(Matriz a, Vector v){

       int n=v.n;  

       Vector b(n);

       for(int i=0; i<n; i++){

           for(int k=0; k<n; k++){

              b.x[i]=b.x[i] + a.x[i][k] * v.x[k];

           }

       }

       return b;

   }



vlinea3d InterseccionTriangulos(poligono3d A, poligono3d B){

//Aquí declaramos las variables necesarias

Matriz X(2);

double rr;

vlinea3d f;

Vector Y(2);



A.plano.coeficientes(A.A, A.B, A.C);   //Calculamos los coeficientes de los planos

B.plano.coeficientes(B.A, B.B, B.C);

A.plano.vector();                      //Calculamos las normales del plano

B.plano.vector();  

   



                     

f.v = CrossProduct(A.plano.v, B.plano.v); //El vector de la solución es este

f.v.normalizar();                         //Se normaliza el vector f.v (pero no el punto f.p)



f.p.x=0;//Y así obtenemos un sistema de dos ecuaciones con dos incógnitas a partir de un sistema de 2 ecuaciones con 3 incógnitas...



/*    Aqui rellenamos (X)*R = Y            */



rr=A.plano.B*B.plano.C+A.plano.C*-B.plano.B;//es el número entre el cual tenemos que dividir cada número de la supuesta inversa de una matriz 2x2 para obtener la matriz inversa real



//La matriz X es la inversa

X.x[0][0] = B.plano.C/rr;

X.x[0][1] = -A.plano.C/rr;

X.x[1][0] = -B.plano.B/rr;

X.x[1][1] = A.plano.B/rr;



Y.x[0]=A.plano.D;

Y.x[1]=B.plano.D;



Y = vector_producto(X,Y);



 f.p.y= Y.x[0];

   f.p.z= Y.x[1];

return f;

}