Archive

Archive for the ‘Matemáticas’ Category

Funciones de Hash

February 20th, 2011 No comments

Una función de hash es aquella que convierte una serie de datos en otra serie de datos más corta (generalmente en una simple variable entera de 32 o 64 bits).
Esto es especialmente útil para poder identificar/representar una secuencia de datos de una forma mucho más rápida y eficiente. Pero hay que tener en cuenta que existen infinitas secuencias de datos que dan como resultado el mismo código de hash por lo que no es un proceso reversible, es decir, es prácticamente imposible obtener la secuencia de datos con la que fue calculado el hash a partir del propio hash.

Por ejemplo, es muy útil para codificar cadenas de caracteres y de esta forma poder realizar búsquedas mucho más rápidamente, o como simple sistema de indexación de datos.
También se utilizan en criptografía y como sistemas de comprobación de integridad de los datos, algunos ejemplos de este tipo de funciones de hash son la CRC, MD5, MD4, Tiger, RMD, SHA, etc.

El principal problema surge cuando es imprescindible que cada secuencia de datos tenga un identificador hash único, si dos o más secuencias tienen el mismo hash se produce una colisión de hash. En este ejemplo las cadenas “Jirafa” y “Tigre” tienen el mismo hash, 02.

Existen infinitas formas de hacer una función de hash y puedes crear la tuya propia fácilmente, el problema es conseguir que la función genere una buena distribución de los hashes y, a ser posible, que produzca las menores colisiones posibles incluso con secuencias de datos muy cortas y similares. Por eso es recomendable utilizar alguna de las muchas funciones de hash que se utilizan desde hace tiempo y están muy probadas.
Aqui dejo un video del MIT de introducción a las técnicas de hashing:

Sin mas historias vamos con algunas de las funciones de hashing de 32 bits más básicas y utilizadas:

Esta surgió de una idea presentada al comité del IEEE POSIX P1003.2 por Glenn Fowler y Phong Vo en 1991, luego Landon Curt Noll la mejoró.
Funciona muy bien, es rápida y produce pocas colisiones. Realiza una buena dispersión de los hashes y va muy bien para hacer hashing en cadenas casi idénticas como URLs, URIs, IPs, nombres de archivo, etc.

1
2
3
4
5
6
7
8
9
10
11
12
UI32 HashFNV(const char *data, UI32 len)
{
  const UI32 fnv_prime = 0x811C9DC5;
  UI32       hash      = 0;
  for(UI32 i=0; i < len; data++, i++)
  {
    hash *= fnv_prime;
    hash ^= (*data);
  }

  return(hash);
}

Esta es del profesor Daniel J. Bernstein, es una de las más eficientes:

1
2
3
4
5
6
7
8
9
10
UI32 HashDJB(const char *data, UI32 len)
{
  UI32 hash = 5381;
  for(UI32 i=0; i < len; data++, i++)
  {
    hash = ((hash << 5) + hash) + (*data);
  }

  return(hash);
}

Esta es de Donald E. Knuth, está en el volumen III de su libro “The Art of Computer Programming”:

1
2
3
4
5
6
7
8
9
10
UI32 HashDEK(const char *data, UI32 len)
{
  UI32 hash = len;
  for(UI32 i=0; i < len; data++, i++)
  {
    hash = ((hash << 5) ^ (hash >> 27)) ^ (*data);
  }
 
  return(hash);
}

Esta está basada en el trabajo de Peter J. Weinberger de los laboratorios AT&T.
En el libro “Compilers (Principles, Techniques and Tools)” se recomienda su uso:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
UI32 HashPJW(const char *data, UI32 len)
{
  const UI32 bits_in_ui32   = (UI32)(sizeof(UI32) * 8);
  const UI32 three_quarters = (UI32)((bits_in_ui32  * 3) / 4);
  const UI32 one_eighth     = (UI32)(bits_in_ui32 / 8);
  const UI32 high_bits      = (UI32)(0xFFFFFFFF) << (bits_in_ui32 - one_eighth);
  UI32       test           = 0;
  UI32       hash           = 0;
  for(UI32 i=0; i < len; data++, i++)
  {
    hash = (hash << one_eighth) + (*data);
    if((test = hash & high_bits) != 0)
    {
      hash = ((hash ^ (test >> three_quarters)) & (~high_bits));
    }
  }

  return(hash);
}

Esta es muy similar a la PJW, pero está optimizada para procesadores de 32 bits.
Se utiliza en la mayoría de los sistemas UNIX:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
UI32 HashELF(const char *data, UI32 len)
{
  UI32 x    = 0;
  UI32 hash = 0;
  for(UI32 i=0; i < len; data++, i++)
  {
    hash = (hash << 4) + (*data);
    if((x = hash & 0xF0000000) != 0)
    {
      hash ^= (x >> 24);
    }
    hash &= ~x;
  }

  return(hash);
}

Esta es del libro “The C Programming Language” de Dennis Ritchie y Brian Kernighan:

1
2
3
4
5
6
7
8
9
10
11
UI32 HashBKDR(const char *data, UI32 len)
{
  const UI32 seed = 131;
  UI32       hash = 0;
  for(UI32 i=0; i < len; data++, i++)
  {
    hash = (hash * seed) + (*data);
  }

  return(hash);
}

Esta es una función de hash del libro “Algorithms in C” de Robert Sedgwicks:

1
2
3
4
5
6
7
8
9
10
11
12
13
UI32 HashRS(const char *data, UI32 len)
{
  UI32 a    = 63689;
  UI32 b    = 378551;
  UI32 hash = 0;
  for(UI32 i=0; i < len; data++, i++)
  {
    hash = hash * a + (*data);
    a    = a * b;
  }

  return(hash);
}

Esta es una función de hash de Justin Sobel:

1
2
3
4
5
6
7
8
9
10
UI32 HashJS(const char *data, UI32 len)
{
  UI32 hash = 1315423911;
  for(UI32 i=0; i < len; data++, i++)
  {
    hash ^= ((hash << 5) + (*data) + (hash >> 2));
  }

  return(hash);
}

Esta realiza una buena distribución de los hash en muchos tipos diferentes de secuencias de datos:

1
2
3
4
5
6
7
8
9
10
UI32 HashSDBM(const char *data, UI32 len)
{
  UI32 hash = 0;
  for(UI32 i=0; i < len; data++, i++)
  {
    hash = (*data) + (hash << 6) + (hash << 16) - hash;
  }

  return(hash);
}

Esta realiza un simple hashing rotativo/aditivo:

1
2
3
4
5
6
7
8
9
10
UI32 HashBP(const char *data, UI32 len)
{
  UI32 hash = 0;
  for(UI32 i=0; i < len; data++, i++)
  {
    hash = (hash << 7) ^ (*data);
  }

  return(hash);
}

Y por último, ésta hecha a partir de las anteriores, también es una mezcla de hash rotativo y aditivo:

1
2
3
4
5
6
7
8
9
10
11
UI32 HashAP(const char *data, UI32 len)
{
  UI32 hash = 0xAAAAAAAA;
  for(UI32 i=0; i < len; data++, i++)
  {
    hash ^= ((i & 1) == 0) ?
              ((hash <<  7) ^ (*data) * (hash >> 3)) :
              (~((hash << 11) + ((*data) ^ (hash >> 5))));
  }

  return(hash);
Categories: C++, Matemáticas Tags: , , ,

Evaluar expresiones matemáticas

July 25th, 2010 No comments

Existen tres notaciones básicas para representar expresiones matemáticas; infija, prefija y postfija.
La infija es la más conocida porque es la que más se utiliza, por ejemplo: (1+2)*(3-4).
En forma prefija se representaría como *+12-34, y en forma postfija sería 12+34-*.

La forma infija resulta muy sencilla para los humanos, estamos acostumbrados a ella y es muy simple. Sin embargo requiere del uso de paréntesis para establecer el orden correcto de las operaciones a realizar, y esto supone un engorro a la hora de programar un evaluador de expresiones.  El “truco” consiste en convertir la expresión de infja a postfija o prefija, de esta forma nos deshacemos de los paréntesis ya que las notaciones prefija y postfija no los necesitan.

Una vez convertida la expresión a notación postfija o prefija es muy fácil crear un árbol binario para calcular el resultado de la expresión.

Los más observadores ya se habrán dado cuenta de que la notación postfija es también conocida como RPN (“Reverse Polish Notation”, notación polaca inversa). Fué inventada por el filósofo y matemático polaco Lukasiewicz precisamente como una notación alternativa para evitar los paréntesis en la expresiones. Más de uno habrá usado esta notación en las calculadoras RPN, especialmente de HP, ¿verdad?

Vamos a ver con un ejemplo muy simple como, usando simplemente un stack,  se puede convertir cualquier expresión infija (por compleja que sea) a notación postfija. Explicaré solamente la conversión de infija a postfija. La conversión a prefija es igual de sencilla pero lo dejo como ejercicio de pasatiempo para el que le apetezca.

Se va cogiendo cada elemento de la expresión (de izquierda a derecha) y se siguen estos pasos:

  1. Si el elemento es un número o variable se copia directamente a la salida.
  2. Si es un paréntesis izquierdo se mete en el stack.
  3. Si es un paréntesis derecho, se van sacando uno a uno todos los elementos del tope del stack (y copiándolos en la salida) hasta encontrar un paréntesis izquierdo en el stack. Finalmente ambos paréntesis se descartan.
  4. Si es un operador y en el tope del stack no hay un operador o hay un operador de menor precedencia entonces se añade el nuevo operador al stack.
  5. Si es un operador y en el tope del stack hay un operador de mayor o igual precedencia entonces se saca el operador del stack y se copia en la salida, y se repite la comparación del operador con el nuevo tope del stack.
  6. Cuando ya no quedan más elementos en la expresión se van sacando uno a uno todos los elementos del stack y se copian en la salida.

Usaré el ejemplo inicial en notación infija: (1+2)*(3-4).

  • ( Se mete en el stack.  
    STACK[ ( ]  SALIDA[ ]

  • 1 Se copia en la salida.
    STACK[ ( ]  SALIDA[ 1 ]

  • + Operador, en el tope del stack no hay otro operador, se mete en el stack.
    STACK[ +( ]  SALIDA[ 1 ]

  • 2 Se copia en la salida.
    STACK[ +( ]  SALIDA[ 12 ]

  • ) Se sacan del stack y se copian en la salida todos los elementos hasta encontrar un paréntesis izquierdo.  
    STACK[ ]  SALIDA[ 12+ ]

  • * Operador, en el tope del stack no hay otro operador, se mete en el stack.  
    STACK[ * ]  SALIDA[ 12+ ]

  • ( Se mete en el stack.  
    STACK[ (* ]  SALIDA[ 12+ ]

  • 3 Se copia en la salida.
    STACK[ (* ]  SALIDA[ 12+3 ]

  • - Operador, en el tope del stack no hay otro operador, se mete en el stack.  
    STACK[ -(* ]  SALIDA[ 12+3 ]

  • 4 Se copia en la salida.  
    STACK[ -(* ]  SALIDA[ 12+34 ]

  • ) Se sacan del stack y se copian en la salida todos los elementos hasta encontrar un paréntesis izquierdo.  
    STACK[ * ]  SALIDA[ 12+34- ]

  • Ya no quedan más elementos en la expresión, se sacan todos del stack y se copian en la salida.  
    STACK[ ]  SALIDA[ 12+34-* ]

El resultado final en notación postfija es 12+34-*

Ahora resulta muy sencillo evaluar la expresión utilizando de nuevo un simple stack.
Vamos cogiendo cada elemento de la expresión en notación postfija (de izquierda a derecha);

  • Si el elemento es un operando (número, variable o constante) se mete en el stack.
  • Si el elemento es un operador (+, -, *, /, sin, sqrt…) se saca del stack tantos elementos como parámetros necesite el operador, se realiza la operación y el resultado se mete en el stack.
  • El proceso se repite hasta que no queden más elementos en la expresión. Al final del proceso habrá en el stack un sólo elemento que corresponde con el resultado final de la evaluación de la expresión.

Veamos como evaluar la expresión 12+34-*

  • 1 Es un operando, se mete en el stack.
    STACK[ 1 ]
  • 2 Es un operando, se mete en el stack.
    STACK[ 21 ]
  • + Es un operador, como + es un operador binario (necesita dos elementos) se sacan del stack y se realiza la operación (1 + 2) y el resultado se mete en el stack.
    STACK[ 3 ]
  • 3 Es un operando, se mete en el stack.
    STACK[ 33 ]
  • 4 Es un operando, se mete en el stack.
    STACK[ 433 ]
  • - Es un operador binario, se sacan dos elementos del stack y se realiza la operación (3 – 4) y el resultado se mete en el stack.
    STACK[ (-1)3 ]
  • * Es un operador binario, se sacan dos elementos del stack y se realiza la operación (3 * (-1)) y el resultado se mete en el stack.
    STACK[ (-3) ]
  • Resultado final: -3

También se puede utilizar un árbol binario para evaluar las expresiones. Para crearlos usamos un stack y vamos cogiendo cada elemento de la expresión en notación postfija (de izquierda a derecha);

  • Si el elemento es un número o variable se crea un nodo con él y se mete en el stack.
  • Si el elemento es un operador binario se crea un nodo con dicho operador, se extraen del stack dos nodos y se le añaden como hijos, finalmente el nuevo nodo se añade al stack.
  • Si el elemento es un operador unuario se crea un nodo con dicho operador, se extrae del stack un nodo y se le añade como hijo, finalmente el nuevo nodo se añade al stack.
  • El proceso se repite hasta que no queden más elementos en la expresión.

Para terminar, simplemente recorriendo recursivamente el árbol en inorden podemos calcular el resultado final de la expresión o incluso reconstruirla en su notación infija:

  1. Imprimimos paréntesis izquierdo.
  2. Visitamos subárbol izquierdo (*).
  3. Imprimimos la raíz.
  4. Visitamos subárbol derecho (*).
  5. Imprimimos paréntesis derecho.

(*): El proceso se repite recursivamente en los pasos 2 y 4.

Categories: C++, Inteligencia Artificial, Matemáticas Tags: , , , , , , , ,

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).

Categories: 3D, C++, Computer Graphics, Física, Matemáticas Tags: , , ,

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.

Categories: 3D, Computer Graphics, Matemáticas Tags: ,

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)))

Categories: 3D, C++, Computer Graphics, Matemáticas, OpenGL Tags: ,

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.

Categories: 3D, C++, Computer Graphics, Matemáticas Tags: , ,