Entradas en "sse"

Intersección Rayo-Triángulo

Hagamos un poco de abstract:

La intersección rayo-triángulo en R3 es un problema clásico que a menudo nos encontramos cuando buscamos soluciones a un sistema que depende de la topología de una malla triangular que interactúa con otra entidad en la escena. Especialmente en aplicaciones y variaciones de raytracing tipo mapeado fotónico; esto es, para renders offline.

its-a-trapframe de  Big Buck Bunny, renderizado con Blender

Hay muchas técnicas para calcularla, una casi omnipresente y eternamente recomendada es la Intersección Rayo-Triángulo Rápida de Mínimo Almacenamiento de Moller y Trumbore (conocida como MT97); basada en el cambio de coordenadas baricéntricas del triángulo en cuestión.

mt97

El paper del último link de MT97 explica bien el algoritmo y proporciona una implementación. Probablemente, es la que necesitáis para vuestro problema.

¿Probablemente?

Lo mejor de MT97 es que nos da una valiosa lección de filosofía. Su defecto, como el de tantas cosas en la vida, es el de su adopción masiva sin planteamiento.

Más que en ningún otro sitio, en este tipo de problemas hay que ser especialmente crítico; sobretodo en las líneas clave que van a suponer más del 70% del código de ejecución (no fuente) del programa.

Así, hace algunos años cuando me encontré con la necesidad de implementar una solución para calcular la intersección rayo/triángulo hice una comparativa (cosas de la intuición) entre MT97 y mi opción alternativa. Sorprendentemente encontré que para mi problema era mucho más rápida (computacionalmente) esta última.

La primera ventaja que intuí antes de hacer la prueba venía por su naturaleza paralelizable. Programándola usando instrucciones SSE y registros XMM* en gran parte del proceso obtenía hasta tres operaciones por instrucción.

El algoritmo es el siguiente:

  1. Calcular la intersección del rayo con el plano que define el triángulo (problema mucho más sencillo y bien conocido -se da en matemáticas de bachillerato <o al menos yo lo dí>-).
  2. Sumar las áreas de los tres triángulos formados al unir los vértices del triángulo original con el punto de intersección (si el rayo es aproximado o totalmente paralelo al plano, terminamos aquí).
  3. Sumar las áreas de los tres triángulos formados y calcular el área del triángulo original.

Se puede demostrar matemáticamente que si la suma de las tres áreas de los formados es igual a la del original, entonces el punto de intersección está dentro del original. Dando como trivial el punto 1, nuestro problema se reduce pues a calcular el área de un triángulo en R3.

La función implementada y con comentarios en español quedaría:

/*
Los registros XMM* son de 128bits, utilizo 4 floats
para llenarlos en una estructura V4.
*/
struct V4
{
   float a;float b;float c;float d;
   V4(float x=0){a=b=c=d=x;}
   V4(float x, float y, float z, float i){a=x;b=y;c=z;d=i;}
};
/* P*: Vértice *(x,y,z) - R: Resultado */
void inline AreaOfATriangle(V4& P1, V4& P2, V4& P3, V4& R)
{
   V4 DIV(0.5);
   //Nota: Las coordenadas .d deben ser 0
   //TODO: Setear a 0 estos últimos 32bits para P1, P2 y P3
      asm volatile
      (
      "MOVUPS %0, %%XMM1\n\t"     //Sube x,y,z de P1 a XMM1
      "MOVUPS %1, %%XMM2\n\t"     //Sube x,y,z de P2 a XMM2
      "MOVUPS %2, %%XMM3\n\t"     //Sube x,y,z de P3 a XMM3
      "MOVUPS %%XMM2, %%XMM4\n\t" //Copia P2 en XMM4
      "SUBPS  %%XMM1, %%XMM4\n\t" //Resta P2-P1 (XMM4)
      "MOVUPS %%XMM3, %%XMM5\n\t" //Copia P3 en XMM5
      "SUBPS  %%XMM1, %%XMM5\n\t" //Resta P3-P1 (XMM5)
      //--- Producto vectorial de P2P1 y P3P1 (determinante 3x3)
      "PSHUFD $0xC9, %%XMM4, %%XMM6\n\t" //P2P1 en XMM6 (Y,Z,X)
      "PSHUFD $0xC9, %%XMM5, %%XMM7\n\t" //P3P1 en XMM7 (Y,Z,X)
      "MULPS  %%XMM4, %%XMM7\n\t" //P2P1 * P3P1(Y,Z,X)
      "MULPS  %%XMM5, %%XMM6\n\t" //P3P1 * P2P1(Y,Z,X)
      "SUBPS  %%XMM7, %%XMM6\n\t" //Resta XMM6 y XMM7
      //--- Módulo del vector resultado
      "MOVUPS %%XMM6, %%XMM7\n\t" //Resultado a XMM7
      "MULPS  %%XMM6, %%XMM7\n\t" //Cuadrado de componentes XMM7
      "HADDPS %%XMM7, %%XMM7\n\t" //Suma parcial de coordenadas
      //... (A+B,C+D,A+B,C+D)
      "HADDPS %%XMM7, %%XMM7\n\t" //Fin la suma
      //... ((A+B)+(C+D),...)
      "SQRTPS %%XMM7, %%XMM7\n\t" //Raíz cuadrada del resultado
      "MOVUPS %3, %%XMM6\n\t"     //Sube 0.5 a XMM6
      "MULPS  %%XMM6, %%XMM7\n\t" //Divide el resultado entre 2
      "MOVUPS %%XMM7, %4"         //Pone el resultado en R {R(0)}
      :"=m"(P1), "=m"(P2), "=m"(P3), "=m"(DIV), "=m"(R)
      :"m"(R)
      );
      return;

Como siempre lo mejor es hacer pruebas en las máquinas que vayan a ejecutar el código. Gran parte de la mejora en velocidad viene del hecho de reducir los accesos a memoria (todo lo que sea pasearse por el bus de la placa base es siempre como subirse a una tortuga…) y no salir de los registros de la CPU.

De esto hace ya algunos años, la extensión SSE4.2 incluye una instrucción para calcular el producto vectorial (nos ahorramos 5 instrucciones en el código dado). Si la integráis con la intersección plano/triángulo y la adaptáis a una CPU contemporánea de gama alta, vuestro raytracer podrá casi volar… 😉

Leer más