Intersecciones entre rayos, esferas y triángulos

La semana pasada estuve teniendo bastantes problemas con la intersección entre rayos y esferas y entre rayos y triángulos. Para no reinventar la rueda opté por usar implementaciones ya existentes, concretamente usé implementaciones de los libros Real-Time Rendering y Real-Time Collision y en principio parecían funcionar, el problema fue cuando añadí transformaciones a los objetos 3D.

Cuando implementé las transformaciones, los objetos no se renderizaban correctamente y estuve bastante tiempo averiguando que es lo que pasaba, pensé que al coger implementaciones de libros tan reconocidos éste no sería el problema, hasta que descartando problemas solo me quedó que los algoritmos de intersección estuviesen mal implementados.

Lo que hice al principio fue buscar por internet otra implementación para no perder mucho el tiempo en implementar yo mismo un algoritmo relativamente sencillo, mi sorpresa fue que tampoco conseguí hacer funcionar las transformaciones en RTal con estos algoritmos. Había probado ya 3 algoritmos distintos y ninguno funcionaba correctamente con las transformaciones, así que finalmente hice lo que debía haber hecho desde un principio, implementarlos yo mismo.

Este artículo va a tratar de como implementar intersecciones entre rayos y esferas y entre rayos y triángulos, intentaré que sea lo más claro posible incluso para aquellas personas que no están muy familiarizadas con la geometría y las matemáticas.

Voy a empezar explicando que es un "rayo" que quizá es un concepto extraño para aquel que no esté familiarizado con el tema.

Un rayo es un concepto geométrico que empieza en un punto y se extiende de manera infinita en una dirección y en un único sentido. Un rayo está muy relacionado con el concepto de linea y de segmento. Una línea se extiende de manera infinita en una dirección y en ambos sentidos, un segmento es un fragmento de recta que une 2 puntos del espacio.

Las definición matemática de rayo, línea y segmento respectivamente es

$R(t) = \textbf{o} + t*\hat{d}$

$L(t) = (1-t)\textbf{a} + t*\textbf{b}$

$S(t) = \textbf{a} + t* (\textbf{b}-\textbf{a})$

En este caso nos interesa solamente el concepto de rayo que es el que se usa en los sistemas de Ray Tracing, el origen del rayo por lo general es la camara o un punto de intersección de la escena y la dirección es la coordenada del pixel que se está calculando, en un futuro artículo hablaré de esto.

Nota: La notación $\textbf{a}$ significa que estamos hablando de un vector y la notación $\hat{a}$ significa que estamos hablando de un vector unitario que se define de la siguiente forma

$\hat{v} = \hat{(x,y,z)} = \frac{\textbf{v}}{||\textbf{v}||}$

donde

$||\textbf{v}|| = \sqrt[2]{x^2 + y^2 + z^2}$

Intersección entre esferas y rayos

Empecemos con las intersecciones entre rayos y esferas.

La definición matemática de una esfera es la siguiente

$(\textbf{p}-\textbf{c})\cdot(\textbf{p}-\textbf{c}) - r^2 = 0$

Donde $\textbf{c}$ es el centro de la esfera y $\textbf{p}$ es un punto cualquiera de su superficie y $r$ es el radio de la esfera, recordemos que el producto punto de dos vectores o dot product es un escalar y se define como

$(a,b,c)\cdot(x,y,z) = a*x + b*y + c*z$

Para encontrar intersecciones entre rayos y esferas basta con definir $\textbf{p}$ de la siguiente forma $\textbf{p} = \textbf{o} + t*\hat{d}$ de esta manera sustituyendo en la ecuación de la esfera $\textbf{p}$ por la ecuación del rayo obtenemos

$(\textbf{o} + t\hat{d} - \textbf{c}) \cdot (\textbf{o} + t\hat{d} - \textbf{c}) - r^2 = 0$

Esta ecuación es una ecuación cuadrática que si simplificamos obtenemos

$t^2(\hat{d}\cdot \hat{d}) + 2t(\hat{d} \cdot (\textbf{o} - \textbf{c})) + (\textbf{o} - \textbf{c}) \cdot (\textbf{o} - \textbf{c}) - r^2 = 0$

Al resorver esta ecuación hay que tener encuenta que las ecuaciones cuadráticas tienen como máximo 2 soluciones (raices) y estas solucciones tienen sentido geométrico.

Si no hay ninguna solución entonces no se ha producido ninguna intersección.
Si hay 2 soluciones significa que el rayo ha entrado en la esfera en un punto del espacio y ha salido en otro punto distinto, en este caso nos interesa el solucion que represente un punto más cercano al origen del rayo.
- Si hay 2 soluciones positivas nos quedamos con la produzca el punto más cercano al origen del rayo.
- Si hay una solución positiva y otra negativa nos quedamos con la solución positiva.
- Si las 2 soluciones son negativas no se produce intersección.
Si solo hay una ecuación significa que el rayo pasa tangente a la esfera.

Una implementación sencilla de las ecuaciones y el proceso explicado sería la siguiente (las 2 primeras líneas del método hit sirven para transformar el rayo en el sistema local de coordenadas de la esfera y así poder realizar de forma sencilla la intersección con esferas transformadas)

Intersección entre rayos y triángulos

Un triángulo está definido por 3 puntos en el espacio llamados vértices (vertex), la unión de estos vértices forman lo que se conoce como aristas (edge) y por lo general un triángulo representa una cara (face).

La intersección entre rayos y triángulos es un poco más compleja pero relativamente sencilla, para saber si un rayo produce una intersección con un triángulo debemos llevar a cabo 2 comprobaciones, la primera es si el rayo intersecta con el plano (infinito) que contiene al triángulo y la segunda comprobar que el punto de intersección (generado por el rayo y el plano) está dentro del triángulo. Así que vamos a dividir esta explicación en 2 partes

Intersección Rayo-Plano

Un plano se define como

$ax + by + cz = d$

También se puede definir de forma vectorial de la siguiente forma

$\hat{n} \cdot (\textbf{x} - \textbf{p}) = 0$

$\hat{n}$ es el vector normal y representa la dirección de la superficie, el producto punto puede representar el coseno del angulo formado por dos vectores, en este caso tiene sentido, si cogemos el vector normal (dirección de la superficie del plano) y un vector del plano, tenemos que el ángulo entre ambos vectores forma $90^\circ$ y $cos(90) = 0$ podemos desarrollar un poco más la fórmula del plano y llegar a

$\hat{n} \cdot \textbf{x} = \hat{n} \cdot \textbf{p}$

Si definimos $\hat{n} \cdot p$ como la distancia desde el origen de coordenadas al plano ( $d$ ) tenemos la siguiente ecuación

$\hat{n} \cdot \textbf{x} = d$

Ya sabemos todo lo que tenemos que saber sobre un plano, ahora para encontrar la intersección entre un rayo y el plano basta con susituir la ecuación del rayo ( $R(t) = \textbf{o} + t*\hat{d}$ ) en la variable $\textbf{x}$ del plano, de esta manera tenemos

$\hat{n} \cdot R(t) = d \longrightarrow \hat{n} \cdot (\textbf{o} + t*\hat{d}) = d \longrightarrow \hat{n} \cdot \textbf{o} + \hat{n}*t \cdot \hat{d} = d$

No confundir $d$ que es la distancia desde el origen de coordenadas al plano con $\hat{d}$ que es la dirección del rayo.

Para resolver esta ecuación tenemos que calcular antes una serie de variables que todavía no conocemos como $\hat{n}$ y $d$ . El vector normal de un plano que contiene un triángulo es el vector normal del triángulo y el vector normal de un triángulo es

$\hat{n} = \frac{(\textbf{b}-\textbf{a}) \times (\textbf{c}-\textbf{a})}{|(\textbf{b}-\textbf{a}) \times (\textbf{c}-\textbf{a})|}$

Donde $\times$ es el producto vectorial y se puede definir de la siguiente forma

$(x,y,z) \times (a,b,c) = \begin{vmatrix} \textbf{i} & \textbf{j} & \textbf{k} \\ a & b & c \\ x & y & z \end{vmatrix} = \textbf{i} \begin{vmatrix} b & c \\ y & z \end{vmatrix} - \textbf{j} \begin{vmatrix} a & c \\ x & z \end{vmatrix} + \textbf{k} \begin{vmatrix} a & b \\ x & y \end{vmatrix}$

Formando $\textbf{i}, \: \textbf{j} \: y \: \textbf{k}$ una base de coordenadas ortonormal y en muchos casos $\textbf{i}=(1,0,0), \: \textbf{j}=(0,1,0) \: y \: \textbf{k}=(0,0,1)$

Teniendo calculado $\hat{n}$ ya solo nos queda calcular la distancia desde el origen $d$ , para ello basta con sustituir todos los datos que tenemos en la ecuación del plano y como punto del plano podemos elegir cualquiera de los vertices del triángulo

$\hat{n} \cdot \textbf{a} = d$

Ahora podemos calcular $t$ para averiguar en que punto el rayo intersecta con el plano

$t = \frac{d-\hat{n} \cdot \textbf{o}}{\hat{n}\cdot \hat{d}}$

¡Ojo! Si $\hat{n}\cdot \hat{d} = 0$ significa que el rayo es paralelo al plano y esto habrá que tenerlo en cuenta a la hora de realizar la implementación.

Ahora podemos calcular el punto donde intersecta el rayo con el plano de la siguiente forma

$\textbf{q} = \textbf{o} + t*\hat{d}$

Comprobando si el punto de intersección está dentro del triángulo

Calculado el punto de intersección (si es que lo hay) ahora tenemos que ver si dicho punto está dentro de todos los bordes del triangulo, si lo está concluiremos que el rayo ha intersectado con el triángulo y procederemos a calcular las coordenadas de intersección.

Para comprobar en lado de un borde se encuentra un determinado punto podemos usar la siguiente fórmula

$[(\textbf{b} - \textbf{a}) \times (\textbf{q} - \textbf{a})] \cdot \hat{n} \geqslant 0$

Deberemos hacer esta comprobación para los 3 bordes

$[(\textbf{b} - \textbf{a}) \times (\textbf{q} - \textbf{a})] \cdot \hat{n} \geqslant 0$

$[(\textbf{c} - \textbf{b}) \times (\textbf{q} - \textbf{b})] \cdot \hat{n} \geqslant 0$

$[(\textbf{a} - \textbf{c}) \times (\textbf{q} - \textbf{c})] \cdot \hat{n} \geqslant 0$

Por último resta saber las coordenadas de la intersección, para ello deberemos usar las coordenadas baricéntricas $\alpha, \: \beta \: y \: \gamma (0 \leqslant \alpha \: \beta \: \gamma \leqslant 1)$

Dichas coordenadas baricéntricas se pueden calcular en relación al area del triángulo

$\alpha = \frac{Area(QBC)}{Area(ABC)} \: \beta = \frac{Area(AQC)}{Area(ABC)} \: \gamma = \frac{Area(ABQ)}{Area(ABC)}$

El producto vectorial tiene un sentido geométrico bastante interesante y es que la magnitud del mismo representa el area formada por los vectores involucrados en el producto vectorial. Teniendo el vector normal $\hat{n}$ podemos expresar la magnitud del producto vectorial de la siguiente forma

$|(\textbf{c}-\textbf{b}) \times (\textbf{q}-\textbf{b})| = [(\textbf{c}-\textbf{b}) \times (\textbf{q}-\textbf{b})] \cdot \hat{n}$

De esta manera tenemos

$\alpha = \frac{[(\textbf{c} - \textbf{b}) \times (\textbf{q}-\textbf{b})] \cdot \hat{n}}{[(\textbf{b} - \textbf{a}) \times (\textbf{c}-\textbf{a})] \cdot \hat{n}} \: \beta = \frac{[(\textbf{a} - \textbf{c}) \times (\textbf{q}-\textbf{c})] \cdot \hat{n}}{[(\textbf{b} - \textbf{a}) \times (\textbf{c}-\textbf{a})] \cdot \hat{n}} \: \gamma = \frac{[(\textbf{b} - \textbf{a}) \times (\textbf{q}-\textbf{a})] \cdot \hat{n}}{[(\textbf{b} - \textbf{a}) \times (\textbf{c}-\textbf{a})] \cdot \hat{n}}$

Con esto hemos hallado las coordenadas baricéntricas, ahora nos falta hallar el punto de intersección del rayo con el triángulo, para ello basta multiplicar las coordenadas baricéntricas con sus respectivos vértices:

$p_{interseccion} = a*\alpha + b*\beta + c*\gamma$

Una implementación de este algoritmo seria la siguiente (de nuevo, las dos primeras líneas del método hit son para transformar el rayo al sistema de coordenadas locales del triángulo, para optimizar se podría hacer la comprobación de la intersección en el sistema de coordenadas global en vez del local)

A continuación una serie de imagenes generadas con RTal que implementa estos algoritmos (https://github.com/RdlP/RTal)