Rellenando triángulos con el Producto Vectorial y su Baricentro

En este post continuaremos la série para desarrollar un motor de renderizado desde 0 utilizando C++. Nuestra misión a lo largo del post será la de rellenar los triángulos generados anteriormente, produciendo superficies 2D sobre las que podremos aplicar colores o texturas.

Primera aproximación realizando un barrido de líneas

Nuestra siguiente tarea será la de rellenar los triángulos dibujados anteriormente, generar superficies 2D a partir de los puntos. Vamos a realizar una aproximación básica que iremos optimizando con el desarrollo del post.

Implementación para dibujar un triángulo

Empezaremos con una implementación simple del triangulo aprovechando la función para dibujar una línea que tenemos del post anterior.

void triangle(Vec2i t0, Vec2i t1, Vec2i t2, TGAImage &image, TGAColor color) { 
    line(t0, t1, image, color); 
    line(t1, t2, image, color); 
    line(t2, t0, image, color); 
}

Que podremos utilizar de la siguiente forma para dibujar triángulos a partir de las aristas (lados):

Vec2i t0[3] = {Vec2i(10, 70),   Vec2i(50, 160),  Vec2i(70, 80)}; 
Vec2i t1[3] = {Vec2i(180, 50),  Vec2i(150, 1),   Vec2i(70, 180)}; 
Vec2i t2[3] = {Vec2i(180, 150), Vec2i(120, 160), Vec2i(130, 180)}; 
triangle(t0[0], t0[1], t0[2], image, red); 
triangle(t1[0], t1[1], t1[2], image, white); 
triangle(t2[0], t2[1], t2[2], image, green);

Obteniendo un resultado como el siguiente:

Implementación para rellenar un triángulo – Superficie 2D

Ahora tenemos que implementar una función que rellene los triángulos. Vamos a realizar un barrido teniendo en cuenta lo siguiente:

• No depende del orden en el que pasemos los vértices o aristas del triángulo.
• Ordenaremos los vertices de los triángulos (con BubbleSort) en función del componente Y.
• Trazaremos una línea horizontal que separe los lados derecho e izquierdo del triángulo.
• Rasterizaremos los píxeles de los lados derecho e izquierdo del triángulo simultáneamente.

void triangle(Vec2i t0, Vec2i t1, Vec2i t2, TGAImage &image, TGAColor color) { 
    // Ordenamos con BubbleSort en función de Y 
    if (t0.y>t1.y) std::swap(t0, t1); 
    if (t0.y>t2.y) std::swap(t0, t2); 
    if (t1.y>t2.y) std::swap(t1, t2); 
    // Marcamos los límites de las dos partes
    line(t0, t1, image, green); 
    line(t1, t2, image, green); 
    line(t2, t0, image, red); 
}

De esta forma, partimos el triángulo en dos partes:

• La parte A (roja) formada por la hipotenusa $t_0$ a $t_2$
• La parte B (verde) formada por los catetos $t_0$ a $t_1$

Y obtenemos un resultado como este

Como veis la parte verde, la parte B, está formada por dos aristas (lados). Vamos a dividir el triangulo con una línea horizontal que pase justo por el vértice común entre esos dos lados. Lo entenderéis mejor con la siguiente imagen que rastericemos.

void triangle(Vec2i t0, Vec2i t1, Vec2i t2, TGAImage &image, TGAColor color) { 
    // Ordenamos los vertices
    if (t0.y>t1.y) std::swap(t0, t1); 
    if (t0.y>t2.y) std::swap(t0, t2); 
    if (t1.y>t2.y) std::swap(t1, t2); 
    int total_height = t2.y-t0.y; 
	// Rasterizamos la mitad del triangulo t0 a t1 en Y
    for (int y=t0.y; y<=t1.y; y++) { 
        int segment_height = t1.y-t0.y+1; 
        float alpha = (float)(y-t0.y)/total_height; 
		// Errores por divison por 0
        float beta  = (float)(y-t0.y)/segment_height; 
        Vec2i A = t0 + (t2-t0)*alpha; 
        Vec2i B = t0 + (t1-t0)*beta; 
        image.set(A.x, y, red); 
        image.set(B.x, y, green); 
    } 
}

En este caso rasterizamos medio triángulo en función de la $Y$ desde $t_0$ a $t_1$. Tendremos algunos problemas con las divisiones por 0. Estos problemas ya los solventamos en el post para dibujar las líneas correctamente. Como veremos en el siguiente paso, no tendremos problema ya que vamos a rellenar las áreas y no nos interesan tanto las lineas de las aristas (lados).

Como veis, hemos partido los triángulos de forma horizontal en dos partes.

Ahora tendremos que rasterizar también el otro lado. Lo haremos simplemente duplicando el loop para los valores de Y entre $t_1$ a $t_2$. Además en este caso dibujaremos todos los pixeles entre ambos lados para ir rellenando el triángulo.

void triangle(Vec2i t0, Vec2i t1, Vec2i t2, TGAImage &image, TGAColor color) { 
    // Ordenamos los vertices en función de Y
    if (t0.y>t1.y) std::swap(t0, t1); 
    if (t0.y>t2.y) std::swap(t0, t2); 
    if (t1.y>t2.y) std::swap(t1, t2); 
    int total_height = t2.y-t0.y; 
	// Mitad inferior del triángulo de t0 a t1
    for (int y=t0.y; y<=t1.y; y++) { 
        int segment_height = t1.y-t0.y+1; 
        float alpha = (float)(y-t0.y)/total_height; 
		// Problema con la división por 0
        float beta  = (float)(y-t0.y)/segment_height; 
        Vec2i A = t0 + (t2-t0)*alpha; 
        Vec2i B = t0 + (t1-t0)*beta; 
		// Rellenamos entre los valores de X de ambos lados  
		// De izquierda a derecha o intercambiamos
        if (A.x>B.x) std::swap(A, B); 
        for (int j=A.x; j<=B.x; j++) { 
            image.set(j, y, color);
        } 
    } 
	// Mitad superior del triángulo de t1 a t2
    for (int y=t1.y; y<=t2.y; y++) { 
        int segment_height =  t2.y-t1.y+1; 
        float alpha = (float)(y-t0.y)/total_height; 
		// Problema con la división por 0
        float beta  = (float)(y-t1.y)/segment_height;
        Vec2i A = t0 + (t2-t0)*alpha; 
        Vec2i B = t1 + (t2-t1)*beta; 
		// Rellenamos entre los valores de X de ambos lados 
		// De izquierda a derecha o intercambiamos
        if (A.x>B.x) std::swap(A, B); 
        for (int j=A.x; j<=B.x; j++) { 
            image.set(j, y, color);
        } 
    } 
}

Ahora nuestros triángulos ya pintan un poco mejor y podemos aplicarles un color. Vamos a simplificar la función anterior juntando los dos bucles en uno para que sea más fácil modificar el código posteriormente.

void triangle(Vec2i t0, Vec2i t1, Vec2i t2, TGAImage &image, TGAColor color) { 
    if (t0.y==t1.y && t0.y==t2.y) return; // Comprobamos que no sea un triangulo degenerado
    // Ordenamos los vertices del triangulo en función de Y
    if (t0.y>t1.y) std::swap(t0, t1); 
    if (t0.y>t2.y) std::swap(t0, t2); 
    if (t1.y>t2.y) std::swap(t1, t2);
	// Rasterizamos el triángulo
    int total_height = t2.y-t0.y; 
    for (int i=0; i<total_height; i++) { 
        bool second_half = i>t1.y-t0.y || t1.y==t0.y; 
        int segment_height = second_half ? t2.y-t1.y : t1.y-t0.y; 
        float alpha = (float)i/total_height; 
        float beta  = (float)(i-(second_half ? t1.y-t0.y : 0))/segment_height;
		// Calculamos los vectores para los dos lados a la vez
        Vec2i A =               t0 + (t2-t0)*alpha; 
        Vec2i B = second_half ? t1 + (t2-t1)*beta : t0 + (t1-t0)*beta; 
		// Rasterizamos la linea de izquierda a derecha respecto X
        if (A.x>B.x) std::swap(A, B); 
        for (int j=A.x; j<=B.x; j++) { 
            image.set(j, t0.y+i, color);
        } 
    } 
}

Sin embargo, si queremos optimizar nuestro Rasterizado de las superficies 2D de los triángulos, esta aproximación solo utiliza 1 hilo de ejecución y podría paralelizarse en 2 ejecutando los dos bucles de la implementación anterior paralelamente.

Sin embargo, tendremos que pensar en una aproximación alternativa que se pueda paralelizar fácilmente en una gran cantidad de núcleos apropiado para las tarjetas gráficas actuales.

Aproximación paralelizable a la rasterización de triángulos

triangulo(vec2 puntos[3]) { 
    vec2 delimitador[2] = encuentra_delimitador(puntos); 
    for (each pixel in delimitador) { 
        if (inside(puntos, pixel)) { 
            put_pixel(pixel); 
        } 
    } 
}

Si ponemos nuestra cabeza en modo paralelizable el pseudo-código anterior parece bastante viable. Si pensamos en una tarjeta gráfica que puede tener cientos o miles de hilos de ejecución paralelos, comprobar si cada pixel de la imagen pertenece o no a una figura es una tarea muy rápida y sencilla.

Coordenadas Baricéntricas para rasterizar superficies 2D de Triángulos

Las Coordenadas Baricéntricas de un triangulo son 3 números reales $\alpha, \beta, \gamma \in [0,1] $ tales que $\alpha + \beta + \gamma = 1$ que permiten parametrizar el interior de un triángulo.

Se considera un triángulo en el plano de vértices $A = (x_A, y_A)$, $B = (x_B, y_B)$ y $C = (x_C, y_C)$. Entonces, cualquier punto del interior del triángulo puede representarse por 3 coordenadas baricéntricas $\alpha, \beta, \gamma$.

La relación entre las coordenadas cartesianas y las baricéntricas es la siguiente:

$$x = \alpha x_A + \beta x_B + \delta x_C$$

$$y = \alpha y_A + \beta y_B + \delta y_C$$

Ahora, volviendo a nuestro caso particular en el que trabajamos con vectores, dado un triangulo $ABC$ y un punto $P$ tendremos que encontrar las coordenadas baricéntricas de $P$ respecto al triángulo $ABC$.

Buscamos tres valores $(1-u-v, u, v)$ de forma que podamos encontrar el punto $P = (1-u-v)A + uB + vC$.

O si lo planteamos desde otro punto de vista, el punto $P$ tendrá las coordenadas $(u, v)$ en la base oblicua $(A, \vec{AB}, \vec{AC})$, de forma que $$P = A + u \vec{AB} + v \vec{AC}$$

Entonces, el problema se reduce a encontrar los valores $u, v$ tales que $$u \vec{AB} + v \vec{AC} + \vec{PA} = 0$$

Ahora tenemos que resolver esta ecuación vectorial, o un sistema de 2 ecuaciones con 2 incógnitas.

$$\left.
u \vec{AB_x} + v \vec{AC_x} + \vec{PA_x} = 0 \atop
u \vec{AB_y} + v \vec{AC_y} + \vec{PA_y} = 0
\right\}$$

Que podemos pasar a su forma matricial:

$$\left.
\begin{bmatrix}
u & v & 1
\end{bmatrix}
\begin{bmatrix}
\vec{AB_x} \\
\vec{AC_x} \\
\vec{PA_x}
\end{bmatrix} = 0
\atop
\begin{bmatrix}
u & v & 1
\end{bmatrix}
\begin{bmatrix}
\vec{AB_y} \\
\vec{AC_y} \\
\vec{PA_y}
\end{bmatrix} = 0
\right\}$$

Por lo tanto estamos buscando un vector $(u, v, 1)$ que sea ortogonal a $(\vec{AB_x}, \vec{AC_x}, \vec{PA_x})$ y  $(\vec{AB_y}, \vec{AC_y}, \vec{PA_y})$.

Utilizando el producto vectorial para calcular las coordenadas baricéntricas de un triángulo

El producto vectorial es una operación entre dos vectores que nos permite obtener un vector en la dirección perpendicular al plano que contiene a dichos vectores. Perfecto para resolver nuestro problema para encontrar las coordenadas baricéntricas del triángulo. Estabamos buscando un vector $(u, v, 1)$ que sea ortogonal a $(\vec{AB_x}, \vec{AC_x}, \vec{PA_x})$ y  $(\vec{AB_y}, \vec{AC_y}, \vec{PA_y})$.

La nueva función de rasterización de la superficie del triángulo:

• Realiza un barrido de todos los pixeles dentro de un cuadro delimitador generado a partir de un triángulo dado.
• Para cada pixel se calculan sus coordenadas baricéntricas:
  • Si tienen algún componente negativo, el píxel se encuentra fuera de la figura.
  • Si tienen ningún componente negativo, el píxel se encuentra dentro de la figura (triángulo) y se debe pintar.

Vec3f barycentric(Vec2i *pts, Vec2i P) { 
	// Producto vectorial
    Vec3f u = cross(Vec3f(pts[2][0]-pts[0][0], pts[1][0]-pts[0][0], pts[0][0]-P[0]), Vec3f(pts[2][1]-pts[0][1], pts[1][1]-pts[0][1], pts[0][1]-P[1])); 
    if (std::abs(u[2])<1) return Vec3f(-1,1,1); // triangulo degenerado devolvemos -1 en vector
    return Vec3f(1.f-(u.x+u.y)/u.z, u.y/u.z, u.x/u.z); 
} 

Con la función anterior calculamos las coordenadas baricéntricas de un punto $P$ respecto al triángulo.

Utilizaremos la siguiente función auxiliar para calcular el producto vectorial:

Vec3f cross(Vec3f a, Vec3f b){
    float x = a.y * b.z - b.y * a.z;
    float y = a.z * b.x - b.z * a.x;
    float z = a.x * b.y - b.x * a.y;
    return Vec3f(x, y, z);
}

Ahora tendremos que calcular el cuadro delimitador alrededor del triángulo. Esto lo haremos con los puntos max y mín equivalentes a la esquina inferior izquierda y superior derecha.

Vec2i bboxmin(image.get_width()-1,  image.get_height()-1); 
Vec2i bboxmax(0, 0); 
Vec2i clamp(image.get_width()-1, image.get_height()-1); 
for (int i=0; i<3; i++) { 
	for (int j=0; j<2; j++) { 
		bboxmin[j] = std::max(0,        std::min(bboxmin[j], pts[i][j])); 
		bboxmax[j] = std::min(clamp[j], std::max(bboxmax[j], pts[i][j])); 
	} 
} 

Iteraremos por las tres esquinas (puntos) del triángulo para calcular el min y max equivalentes a la esquina inferior izquierda y superior derecha del cuadro delimitador.

Solo nos queda barrer el cuadro delimitador que hemos generado y comprobar cada uno de sus puntos en coordenadas baricéntricas del triángulo para ver si pertenecen o no.

Vec2i P; 
for (P.x=bboxmin.x; P.x<=bboxmax.x; P.x++) { 
	for (P.y=bboxmin.y; P.y<=bboxmax.y; P.y++) { 
		Vec3f bc_screen  = barycentric(pts, P); 
		if (bc_screen.x<0 || bc_screen.y<0 || bc_screen.z<0) continue; 
		image.set(P.x, P.y, color); 
	} 
} 

Y generaremos de nuevo un triángulo con está implementación:

Vec2i pts[3] = {Vec2i(10,10), Vec2i(100, 30), Vec2i(190, 160)}; 
triangle(pts, frame, TGAColor(255, 0, 0, 0)); 

Ahora que ya sabemos como rasterizar un triángulo de forma paralelizable vamos a aplicarlo al modelo 3D que renderizamos en el post anterior.

Renderizado de un modelo 3D con triángulos de colores – Flat Shading

En el post anterior vimos como dibujar los triángulos que forman un modelo 3D. Ahora que hemos aprendido a pintarlos por dentro de forma paralelizable vamos a intentar juntar ambas implementaciones.

1. Adaptamos los puntos del modelo a la escala de la imagen que vamos a renderizar.
2. Pintamos cada uno de los triángulos con un color aleatorio.

Se puede entender mejor con el siguiente código que utilizaremos para rasterizar todos los triángulos del modelo, uno a uno:

for (int i=0; i<model->nfaces(); i++) { 
    std::vector<int> face = model->face(i); 
    Vec2i screen_coords[3]; 
    for (int j=0; j<3; j++) { 
        Vec3f world_coords = model->vert(face[j]); 
        screen_coords[j] = Vec2i((world_coords.x+1.)*width/2., (world_coords.y+1.)*height/2.); 
    } 
    triangle(screen_coords, image, TGAColor(rand()%255, rand()%255, rand()%255, 255)); 
}

Obteniendo una imagen como la siguiente:

Como vemos, el resultado es algo extraño. Nos falta un elemento fundamental que es la iluminación y las sombras.

Vamos a quitar todos esos colores aleatorios que habíamos generado y aplicaremos una iluminación y sombreado.

Aproximación a un renderizado con iluminación frontal

En este caso nos vamos a apoyar en un par de reglas sobre la iluminación:

• Con un foco lumínico que emite la misma intensidad lumínica, el polígono se ilumina más por los lados que son ortogonales a la dirección de la luz.
• La cara del polígono tendrá 0 iluminación si es totalmente paralela a la dirección de la luz.

La intensidad lumínica que recibe una cara de un polígono es igual al producto escalar del vector iluminación y la normal de la cara de un polígono. En nuestro caso, calcularemos la normal de un triángulo, que es el producto vectorial de sus lados.

En este caso trabajaremos con una rasterización lineal de los colores. Sin embargo, el color (128, 128, 128) no tiene ni la mitad de brillo que (255, 255, 255). Ignoraremos la corrección de gamma y trabajaremos con colores con un brillo incorrecto.

Calcularemos la normal de cada triángulo y la utilizaremos para calcular la intensidad de luz que nos servirá para colorear cada triángulo. Lo entenderemos mas con la implementación:

Vec3f light_dir(0,0,-1);
    for (int i=0; i<model->nfaces(); i++) {
        std::vector<int> face = model->face(i);
        Vec2i screen_coords[3];
        Vec3f world_coords[3];
        for (int j=0; j<3; j++) {
            Vec3f v = model->vert(face[j]);
            screen_coords[j] = Vec2i((v.x+1.)*width/2., (v.y+1.)*height/2.);
            world_coords[j] = v;
        }
        Vec3f n = (world_coords[2]-world_coords[0]) ^ (world_coords[1]-world_coords[0]);
        n.normalize();
        float intensity = n*light_dir;
        if (intensity>0) {
            triangle(screen_coords, image, TGAColor(intensity*255, intensity*255, intensity*255, 255));
        }
    }

En la clase Vec3 añadiremos la implementación del producto vectorial y la normalización:

inline Vec3<t> operator ^(const Vec3<t> &v) const { 
	return Vec3<t>(y*v.z-z*v.y, z*v.x-x*v.z, x*v.y-y*v.x); 
}
float norm () const { 
	return std::sqrt(x*x+y*y+z*z); 
}
Vec3<t> & normalize(t l=1) { 
	*this = (*this)*(l/norm());
	return *this; 
}

Obteniendo un resultado final como este. Como veis, en la boca y los ojos obtenemos un rasterizado «raro». Esto lo mejoraremos en el siguiente post en el que utilizaremos un buffer-z. En este caso, para acelerar el renderizado, aquellas caras que  estén «de espaldas» a la fuente de luz no las rasterizamos. Hacemos esto comprobando si el producto es positivo. Esta técnica se denomina Back-face culling.

Como siempre podeis encontrar la fuente de este código en el GitHub del autor.