Archive

Archive for the ‘3D’ Category

Interpolando que es gerundio

June 20th, 2010 No comments

Hay pocas cosas tan útiles en programación como la interpolación, de hecho son tantos los casos en los que resulta útil (o incluso esencial) que ni me voy a molestar en mencionarlos.
Directamente dejo aqui un par de ejemplos, el primero es una función de interpolación lineal:

1
2
3
4
float InterpLineal(float tiempo, float inicio, float fin)
{
  return((1.0f - tiempo) * inicio + tiempo * fin);
}

Esta es una función de interpolación cuadrática (“Easy In”):

1
2
3
4
float InterpCuadIn(float tiempo, float inicio, float fin)
{
  return(InterpLineal(tiempo * tiempo, incio, fin));
}

Y por último una interpolación cuadrática (“Easy In-Out”):

1
2
3
4
5
6
7
8
9
10
11
12
float InterpCuadInOut(float tiempo, float inicio, float fin)
{
  float medio = (inicio + fin) * 0.5f;
  tiempo *= 2.0f;
  if(tiempo <= 1.0f)
  {
    return(InterpLineal(tiempo * tiempo, inicio, medio));
  }
  tiempo -= 1.0f;

  return(InterpLineal(tiempo * tiempo, medio, fin));
}

En los tres casos se fija el rango a interpolar con los parámetros inicio y fin, y con el parámetro tiempo se indica el punto de interpolación (por supuesto se asume que tiempo tendrá un valor comprendido entre 0.0 y 1.0).

Transformaciones inversas, inversión de matrices

August 17th, 2009 No comments

En programación 3D en tiempo real es imprescindible optimizar todo lo posible, y la inversión de matrices homogéneas es uno de los casos complejos (aparentemente) en donde se puede conseguir una excelente optimización siguiendo las reglas del álgebra matricial y conociendo cómo ha sido generada la matriz a invertir (conociendo su estructura).

Nota: En este post usaré el formato fila (pre-multiplicación) para representar las matrices y la formulación.

Vamos a ver primero la inversión aislada de los tres casos básicos, matriz de rotación, escala y traslación:

Si la matriz de rotación R es ortogonal (si sus tres ejes son perpendiculares entre si), entonces su inversa es simplemente su traspuesta:

R=\begin{bmatrix}RXx&RXy&RXz&0\\RYx&RYy&RYz&0\\RZx&RZy&RZz&0\\0&0&0&1\end{bmatrix}

R^{-1}=R^{T}=\begin{bmatrix}RXx&RYx&RZx&0\\RXy&RYy&RZy&0\\RXz&RYz&RZz&0\\0&0&0&1\end{bmatrix}

La inversa de una matriz de escala S es muy sencilla:

S=\begin{bmatrix}Sx&0&0&0\\0&Sy&0&0\\0&0&Sz&0\\0&0&0&1\end{bmatrix}

S^{-1}=\begin{bmatrix}\frac{1}{Sx}&0&0&0\\0&\frac{1}{Sy}&0&0\\0&0&\frac{1}{Sz}&0\\0&0&0&1\end{bmatrix}

Y la inversa de una matriz de traslación T también es muy simple:

T=\begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&1&0\\Tx&Ty&Tz&1\end{bmatrix}

T^{-1}=\begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&1&0\\-Tx&-Ty&-Tz&1\end{bmatrix}

Ahora veamos un caso más práctico, cómo invertir una matriz de transformación RT (es decir, una matriz que primero realiza una rotación y luego una traslación). Es imprescindible saber siempre el orden de las transformaciones contenidas en una matriz para poder optimizar el cálculo de su inversa ya que vamos a optimizar los cálculos por partes y siguiendo una regla muy simple; la inversa del producto de dos matrices es igual al producto invertido de las inversas de dichas matrices: (A*B)^{-1}=B^{-1}*A^{-1}

Con esta regla tan simple podemos invertir una matriz RT multiplicando la inversa de la traslación por la inversa de la rotación: (R*T)^{-1}=T^{-1}*R^{-1} y como ya sabemos que la inversa de una matriz ortogonal es su traspuesta la fórmula final será: (R*T)^{-1}=T^{-1}*R^{T}

Con lápiz y papel calculamos la estructura final de la matriz inversa de una transformación RT:

\large(R*T)^{-1}=\begin{bmatrix}RXx&RYx&RZx&0\\RXy&RYy&RZy&0\\RXz&RYz&RZz&0\\-(RXx*Tx+RXy*Ty+RXz*Tz)&-(RYx*Tx+RYy*Ty+RYz*Tz)&-(RZx*Tx+RZy*Ty+RZz*Tz)&1\end{bmatrix}

Para invertir una matriz TR es incluso más simple: (T*R)^{-1}=R^{T}*T^{-1}

(T*R)^{-1}=\begin{bmatrix}RXx&RYx&RZx&0\\RXy&RYy&RZy&0\\RXz&RYz&RZz&0\\-Tx&-Ty&-Tz&1\end{bmatrix}

Ahora veremos cómo invertir matrices que incluyen transformaciones de escala isométrica (que significa que se ha aplicado la misma escala a los tres ejes), empecemos por la clásica transformación SRT (donde primero se escala, luego se rota y finalmente se traslada). Como siempre su inversa será: (S*R*T)^{-1}=T^{-1}*R^{T}*S^{-1}

Volvemos a coger el lápiz y calculamos:

\large(S*R*T)^{-1}=\begin{bmatrix}\frac{RXx}{S}&\frac{RYx}{S}&\frac{RZx}{S}&0\\\frac{RXy}{S}&\frac{RYy}{S}&\frac{RZy}{S}&0\\\frac{RXz}{S}&\frac{RYz}{S}&\frac{RZz}{S}&0\\-\frac{RXx*Tx+RXy*Ty+RXz*Tz}{S}&-\frac{RYx*Tx+RYy*Ty+RYz*Tz}{S}&-\frac{RZx*Tx+RZy*Ty+RZz*Tz}{S}&1\end{bmatrix}

La misma inversa SRT pero de una matriz con escalado anisométrico (escalado no idéntico en los tres ejes) sería:

\large(S*R*T)^{-1}=\begin{bmatrix}\frac{RXx}{Sx}&\frac{RYx}{Sy}&\frac{RZx}{Sz}&0\\\frac{RXy}{Sx}&\frac{RYy}{Sy}&\frac{RZy}{Sz}&0\\\frac{RXz}{Sx}&\frac{RYz}{Sy}&\frac{RZz}{Sz}&0\\-\frac{RXx*Tx+RXy*Ty+RXz*Tz}{Sx}&-\frac{RYx*Tx+RYy*Ty+RYz*Tz}{Sy}&-\frac{RZx*Tx+RZy*Ty+RZz*Tz}{Sz}&1\end{bmatrix}

Para terminar veremos la inversa de una matriz que represente una transformación RTS. La inversa con escalado isométrico es:

\large(R*T*S)^{-1}=\begin{bmatrix}\frac{RXx}{S}&\frac{RYx}{S}&\frac{RZx}{S}&0\\\frac{RXy}{S}&\frac{RYy}{S}&\frac{RZy}{S}&0\\\frac{RXz}{S}&\frac{RYz}{S}&\frac{RZz}{S}&0\\-(RXx*Tx+RXy*Ty+RXz*Tz)&-(RYx*Tx+RYy*Ty+RYz*Tz)&-(RZx*Tx+RZy*Ty+RZz*Tz)&1\end{bmatrix}

Y la inversa de una matriz RTS con escala anisométrica es:

\large(R*T*S)^{-1}=\begin{bmatrix}\frac{RXx}{Sx}&\frac{RYx}{Sy}&\frac{RZx}{Sz}&0\\\frac{RXy}{Sx}&\frac{RYy}{Sy}&\frac{RZy}{Sz}&0\\\frac{RXz}{Sx}&\frac{RYz}{Sy}&\frac{RZz}{Sz}&0\\-(RXx*Tx+RXy*Ty+RXz*Tz)&-(RYx*Tx+RYy*Ty+RYz*Tz)&-(RZx*Tx+RZy*Ty+RZz*Tz)&1\end{bmatrix}

Si teneis dificultades para ver alguna de las matrices pulsad con el botón derecho del ratón sobre ella y elegid Abrir/Ver imagen.

Formato interno de las matrices en OpenGL

August 15th, 2009 No comments

Desde los inicios de OpenGL existe una confusión bastante importante sobre cómo administra internamente las matrices. Unos piensan que las almacena en formato fila y otros piensan que lo hace en formato columna. Pues lo cierto es que ni una cosa ni la otra!
OpenGL almacena las matrices en un array unidimensional de 16 elementos y es el programador quien decide cómo prefiere interpretarlas.

Quedará más claro con un ejemplo:

1
2
GLfloat m[16];
glGetFloatv(GL_MODELVIEW_MATRIX, m);


Con el código anterior leemos la matriz actual MODELVIEW y la almacenamos en el array m.

¿Y cómo están estructurados los elementos de la matriz? Pues muy sencillo, secuencialmente:

m = \{m0,m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15\}

Ahora el programador puede decidir interpretarlo como una matriz de tipo columna:

\begin{bmatrix}m0&m4&m8&m12\\m1&m5&m9&m13\\m2&m6&m10&m14\\m3&m7&m11&m15\end{bmatrix}

O como una matriz de tipo fila:

\begin{bmatrix}m0&m1&m2&m3\\m4&m5&m6&m7\\m8&m9&m10&m11\\m12&m13&m14&m15\end{bmatrix}

Matemáticamente es indiferente cómo prefieras interpretarlas, pero recuerda que si utilizas el formato columna debes realizar post-multiplicaciones, y si usas el formato fila pre-multiplicaciones.

Supongamos que tienes una matriz R que representa una rotación y la matriz T que realiza una traslación, y quieres combinarlas en una matriz C, de forma que la matriz de trasformación resultante realice primero la rotación y luego la traslación. Si usas el formato fila tendrás que hacer una pre-multiplicación: C=R*T (Esto se lee: la matriz R pre-multiplica a la matriz T).
Pero si usas el formato columna el orden se invierte ya que es una post-multiplicación: C=T*R (Y esto se lee: la matriz R post-multiplica a la matriz T).

Otro ejemplo; quieres transformar el vector V=\{x,y,z,0\} por la matriz M para obtener el vector transformado V'.
Si la matriz M es de tipo fila pre-multiplicas: V'=V*M
En cambio, si es de tipo columna post-multiplicas: V'=M*V

En algebra matricial siempre se multiplican filas por columnas, por lo tanto si la matriz M es de tipo fila V y V' deben ser interpretados como un vectores fila, y si M es de tipo columna entonces V y V' deberán ser interpretados como vectores columna.

Un último ejemplo; transforma el vector V por las matrices A,B,C (en éste orden).
La pre-multiplicación sigue el orden natural de las transformaciones por lo que resulta más intuitiva: V'=V*A*B*C
La post-multiplicación sería el proceso inverso: V'=C*B*A*V

Para terminar un truco para no liarte; cuando utilices matrices y vectores columna (post-multiplicación) debes interpretar las transformaciones como series de funciones: V'=C(B(A(V)))

Raíz cuadrada, optimización redonda

May 7th, 2009 No comments

En programación 3D hay pocas funciones tan necesarias como la raíz cuadrada. ¿Quién no recuerda el famoso teorema de Pitágoras? El cuadrado de la hipotenusa es igual a la suma de los cuadrados de los catetos h^{2}=x^{2} + y^{2} por lo tanto la hipotenusa será h=\sqrt{x^{2} + y^{2}}

Esto lo podemos emplear con cualquier número de dimensiones, en 3D tendremos tres coordenas \left\{x,y,z\right\} que definen un punto o vector en un espacio tridimensional, aplicando el mismo teorema podemos saber el módulo (longitud) del vector m=\sqrt{x^{2}+y^{2}+z^{2}} y a partir del módulo podemos, por ejemplo, normalizarlo (hacer que el vector sea de módulo 1).

Un ejemplo práctico:

1
2
3
4
5
6
7
8
9
10
11
12
V3 &NormalizarVector(V3 &v)
{
  float modulo = sqrt( v.x*v.x + v.y*v.y + v.z*v.z);

  assert(modulo != 0.0f);

  v.x /= modulo;
  v.y /= modulo;
  v.z /= modulo;

  return(v);
}

La función anterior para normalizar vectores es perfectamente correcta pero se puede optimizar bastante. En primer lugar calcular raíces cuadradas tiene un coste computacional bastante importante y en aplicaciones 3D como juegos se usa constantemente por lo que es interesante optimizarla todo lo posible. Además de la raíz cuadrada también podemos sustituir las tres divisiones por tres multiplicaciones, ¿cómo lo hacemos? usando la inversa de la raíz cuadrada \frac{1}{\sqrt{x}} que nos permitirá multiplicar en lugar de dividir ya que la división está implícita en el cálculo.

¿Y cómo optimizamos la raíz cuadrada? Existen muchas formas de hacerlo pero vamos a ver la más rápida y sencilla, la raíz cuadrada inversa rápida (Fast InvSqrt):

1
2
3
4
5
6
7
8
9
10
float InvSqrt (float x)
{
  float x_div2 = x * 0.5f;         // Dividimos x entre 2
  int i = *(int*) &x;              // Pasamos de float a int (ambos de 32 bits)
  i = 0x5f3759df - (i >> 1);       // Hacemos una primera aproximación del resultado
  x = *(float*) &i;                // Pasamos de int a float
  x = x * (1.5f - x_div2 * x * x); // Recalculamos para aproximar más el resultado

  return(x);
}

Esta función se aprovecha de la forma interna de trabajo de la unidad de punto flotante (normativa IEEE 754)  para calcular una buena aproximación del resultado exacto. No es más que una simple aproximación por el método Newton.
Este código se hizo muy popular por haber sido utilizado en el Quake 3 y aunque su origen no está nada claro parece que surgió de SGI (Silicon Graphics) hace bastantes años.

La nueva función de normalización vectorial quedaría de la siguiente manera:

1
2
3
4
5
6
7
8
9
10
V3 &NormalizarVector(V3 &v)
{
  float m = InvSqrt( v.x*v.x + v.y*v.y + v.z*v.z);

  v.x *= m;
  v.y *= m;
  v.z *= m;

  return(v);
}

El incremento de velocidad conseguido es muy importante pero también hay que tener en cuenta que el resultado es una aproximación y que este código sólo puede usarse en procesadores cuyas unidades de punto flotante sigan la normativa IEEE 754.