lunes, 22 de noviembre de 2010

GPU-RT: Geometria constructiva

La geometria constructiva (GC) es un metodo para crear objetos geometricos mas complejos a partir de geometrias basicas mas simples, (usualmente esferas, conos, cajas, etc) aplicando operaciones de interseccion, union y diferencia.
Por ejemplo la interseccion de un cono con un plano puede generar una parabola.

Para un modelador de mallas que se base en triangulos, no es nada trivial generar una malla a partir de 2 mallas originales y una operacion de interseccion o diferencia. Particularmente calcular la frontera entre los 2 objetos originales que se intersectan lleva a varias complicaciones, y la precision siempre es un problema. (*)

Afortundamente para un ray-tracer la GC es una operacion trivial, (o relativamente simple).

El caso de la interseccion A ^ B, si el rayo intersecta al objeto A en el pto Ip, solo hay que chequear que dicho pto pertenezca al volumen de B, y viceversa, si el rayo intersecta en el B, hay que chequear que el pto este en el interior de A.

Ejemplo para la interseccion con 3 esferas:
// Geometria constructiva
// Objeto Interseccion M Interseccion M2
float4 M = tex2D( g_samObj2, float2(j+0.5,i+0.5)/OBJ_SIZE);
float4 M2 = tex2D( g_samObj3, float2(j+0.5,i+0.5)/OBJ_SIZE);
if(fl)
{
float3 pt = LF+D*t0;
float dist = distance(pt,float3(M.x,M.y,M.z));
// interseccion
if(dist>M.w)
fl = false;
else
{
dist = distance(pt,float3(M2.x,M2.y,M2.z));
if(dist>M2.w)
fl = false;
}
}
if(fl2)
{
float3 pt = LF+D*t1;
float dist = distance(pt,float3(M.x,M.y,M.z));
// interseccion
if(dist>M.w)
fl2 = false;
else
{
dist = distance(pt,float3(M2.x,M2.y,M2.z));
if(dist>M2.w)
fl2 = false;
}
}







El caso de la diferencia si A y B son 2 objetos el objeto A-B, se puede resolver en el tracing de la siguiente forma.
Supongamos que el rayo intersecto con el objeto A en el pto Ip, solo hay que chequear que Ip este incluido en el volumen de A. Si no llega a ser el caso, se anula esa interseccion rayo-A.
Usualmente A-B se completa agregando tambien B ^ A, para que dibuje el objeto entero.




Implementar la GC en forma general requiere hallar una formula o un metodo de determinar si un punto cualquiera esta o no dentro del volumen del objeto.
Para el caso de la esfera eso es trivial, y en general en el mismo momento que uno analiza las prop. geometricas de un cuerpo para hallar la formula de la interseccion con el rayo, con un poco mas de trabajo se puede llegar al calculo de inclusion del punto.
Si la formula analitica es muy dificil de hallar, siempre se puede calcular de la siguiente forma, si el cuerpo es "cerrado" e decir no tiene ningun agujero:
desde el pto en cuestion se traza un rayo hacia cualquier direccion, y se cuenta la cantidad de veces que intersecta sobre si mismo. Si esa cantidad es impar el punto estaba adentro del cuerpo. Justamente, la interseccion entre el rayo y el cuerpo ya la resuelta para el ray-tracing. (De lo contrario ese cuerpo no seria una primitiva de nuestro tracer)
Con lo cual implementar la GC en un ray-tracer es bastante directo.


(*) Hace unos años estaba implementando la Geometría Constructiva en el motor grafico que utilizo para mis sistemas, que esta basado en mallas. Internamente todos los objetos se convierten en una malla de triangulos, que luego se envian ahora a la GPU, antes al GDI del windows.
Generar la malla resultante fue sumamente dificil, y si bien a nivel grafico parecia funcionar bien, (siempre se puede ajustar el nivel de precision para que quede bien), a nivel analitico la frontera no quedaba perfectamente cerrada. Como la malla de salida se tenia que enviar en formato STL a una impresora 3d, el driver de la impresora no aceptaba huecos en la figura, asi como tapoco triangulos sueltos.
Programas como el Rhino que trabajan sobre Nurbs tienen implementado en su nucleo todas las operaciones de GC directamente sobre las curvas, y generan las mallas solo en el momento que se necesitan y con la precision requerida.

Ejemplo de un Toro menos otro toro, donde se ve que graficamente esta bien:



Sin embargo, de cerca la frontera, tiene huecos en algunos lugares. Ademas se ve como hay que empezar a tessellar los triangulos en la interseccion entre ambas figuras:

domingo, 21 de noviembre de 2010

GPU-RT: Estabilidad numerica

Siguiendo con la idea de incorporar otras primitivas al ray-tracer, una muy interesante es el toro (torus), o dona. Esta superficie es muy estudiada en topologia, y se obtiene a partir de un rectangulo uniendo primero los 2 lados horizontales formando un cilindro, y luego los 2 verticales formando una dona.
Las ecuaciones parametricas del toro de radio mayor R y radio menor r son:

x(u, v) = (R + r cos v) cos u
y(u, v) = (R + r cos v) sen u
z(u, v) = r sen v

Para una descripcion mas detallada del toru sus ecuaciones y propiedades:
http://es.wikipedia.org/wiki/Toro_(geometr%C3%ADa)

La interseccion entre el rayo y el toro requiere resolver la raiz de un polinomio de grado 4:

float AA = LF.x - toro.x;
float BB = LF.y - toro.y;
float CC = LF.z - toro.z;
float C = AA*AA+BB*BB+CC*CC-(toro_rx*toro_rx + toro_ry*toro_ry);
float B = 2*(D.x*AA+D.y*BB+D.z*CC);
float a2 = toro_rx*toro_rx;
float b2 = toro_ry*toro_ry;

// coefientes ecuacion cuarta
float r0 = 2*B;
float r1 = B*B + 2*C + 4*a2*D.z*D.z;
float r2 = 2*B*C + 8*a2*CC*D.z;
float r3 = C*C+ 4*a2*(CC*CC - b2);

entonces, la ecuacion que hay que resolver es
x^4 + r0 x^3 + r1 x^2 + r2 x + r3 = 0

Si metemos esto en un PS, nos queda asi:

// ecuaciones cuarticas
double solve_cubic(double a, double b, double c);

int solve_quartic(double a, double b, double c, double d, inout float results[4])
{
int num_roots = 0;
double a2 = a*a;
double k = a2 - 4*b;
double twice_y0 = solve_cubic(-b, a*c - 4*d,-k*d - c*c);
double alpha2 = k/4 + twice_y0;
if (alpha2 >=-EPSILON)
{
double alpha,beta;
double y0 = twice_y0/2;
if (abs(alpha2)<EPSILON)
{
alpha = 0.0;
beta = sqrt(y0*y0 - d);
}
else
{
alpha = sqrt(alpha2);
beta = (a*y0 - c)/(2*alpha);
}

double BB = a/2 - alpha;
double CC = y0 - beta;
double disc_sqr = BB*BB - 4*CC;
double disc;
if (abs(disc_sqr)<EPSILON)
{
results[0] = -BB / 2;
num_roots++;
}
else if (disc_sqr > 0)
{
disc = sqrt(disc_sqr);
results[0] = (-BB + disc)/2;
results[1] = (-BB - disc)/2;
num_roots+=2;
}

BB = a/2 + alpha;
CC = y0 + beta;
disc_sqr = BB*BB - 4*CC;
if (abs(disc_sqr)<EPSILON)
{
if(num_roots==0)
results[0] = -BB / 2;
else
if(num_roots==1)
results[1] = -BB / 2;
else
results[2] = -BB / 2;
num_roots++;
}
else if (disc_sqr > 0)
{
disc = sqrt(disc_sqr);
if(num_roots==0)
{
results[0] = (-BB + disc)/2;
results[1] = (-BB - disc)/2;
}
else
if(num_roots==1)
{
results[1] = (-BB + disc)/2;
results[2] = (-BB - disc)/2;
}
else
{
results[2] = (-BB + disc)/2;
results[3] = (-BB - disc)/2;
}
num_roots+=2;
}
}

return num_roots;
}


double solve_cubic(double a, double b, double c)
{
double Q = (a*a - 3*b)/9;
double R = (2*a*a*a - 9*a*b + 27*c)/54;
double Q3 = Q*Q*Q;
double disc = R*R - Q3;
double k2 = -a/3;
if (disc < -EPSILON)
{
// 3 raices reales
double theta = acos(R/sqrt(Q3));
double k1 = -2*sqrt(Q);
return k1*cos(theta/3) + k2;
//results[1] = k1*cos((theta + 2*M_PI)/3) + k2;
//results[2] = k1*cos((theta + 4*M_PI)/3) + k2;
//return 3;
}
else
{
// una raiz real
double disc2 = sqrt(disc) + abs(R);
double cuberoot = pow(disc2,1/3.0);
if (R < 0)
return (cuberoot + Q/cuberoot) + k2;
else
return -(cuberoot + Q/cuberoot) + k2;
}
}


void PS_RayTracing2( float4 Diffuse:COLOR0, float2 Tex : TEXCOORD0,
out float4 Color : COLOR0 )
{
// Calculo la direccion del rayo
float3 D = normalize(g_vViewDir + g_vDx*(2*Tex.x-1) + g_vDy*(1-2*Tex.y));
float R = 100000;
float3 Ip = 0;

float AA = LF.x - toro.x;
float BB = LF.y - toro.y;
float CC = LF.z - toro.z;
float C = AA*AA+BB*BB+CC*CC-(toro_rx*toro_rx + toro_ry*toro_ry);
float B = 2*(D.x*AA+D.y*BB+D.z*CC);
float a2 = toro_rx*toro_rx;
float b2 = toro_ry*toro_ry;

// coefientes ecuacion cuartica
float r0 = 2*B;
float r1 = B*B + 2*C + 4*a2*D.z*D.z;
float r2 = 2*B*C + 8*a2*CC*D.z;
float r3 = C*C+ 4*a2*(CC*CC - b2);

// Resuelvo la ecuacion cuartica
float roots[] = {0,0,0,0};
int cant = solve_quartic(r0,r1,r2,r3,roots);

if(cant>0)
Color.rgb = LF+D*roots[0]; // punto de interseccion
else
Color.rgb = 0;
Color.a = 1;
}




Desafortundamente, al correr este shader, comienzan a aparecer "huecos" por algunos lugares, en ciertos puntos donde por errores de redondeo o falta de precision, el algoritmo se vuelve numericamente inestable:





El algoritmo se hace mucho mas inestable en la medida que el toro se aleja o se alinea en ciertos angulos:



Como la solucion analitica sufre de problemas de estabilidad, probe de implementar una solucion numerica. El metodo de Laguerre para encontrar raices de polinomios tiene la ventaja de ser muy facil de implementar, y no requiere elegir un punto inicial especial (como el metodo dew Newton), con lo cual siempre converge (Se llaman metodos con convergencia Global, para distinguirlos del metodo de Newton que solo goza de convergencia Local)

http://en.wikipedia.org/wiki/Laguerre's_method

La implementacion en el PS:


double Laguerre(double x,double a, double b, double c, double d)
{
for(int k=0;k<20;++k)
{
double x2 = x*x;
double x3 = x2*x;
double x4 = x3*x;

double p0 = x4+a*x3+b*x2+c*x+d; // p(x)

if(abs(p0)<EPSILON) // ya encontre un raiz con suf. precision
k = 21;
else
{
double p1 = 4*x3+3*a*x2+2*b*x+c; // p'(x)
double p2 = 12*x2+6*a*x+2*b; // p''(x)
double G = p1/p0;
double H = G*G - p2/p0;
double disc = sqrt(abs(3*(4*H-G*G)));
double a = G>0?4/(G+disc):4/(G-disc);
x -= a;
}
}
return k>20?x:0;
}



Al correr el PS, se obtienen tambien puntos negros donde el rayo pierde la superficie del toro, con lo cual tampoco resulta numericamente estable:



Otra vez, la cantidad de pasos (ahora puesta en 20) es arbitraria, pero en este caso es muy dificil de encontrar un valor adecuado para todos los casos y que no sea demasiado alto para que no se haga tan lento.

Re-escribiendo las ecuaciones del toro, para que sea una funcion de 2 variables:
....
float d = R - sqrt(x*x + y*y);
float disc = r*r-d*d;
return sqrt(disc);
....


Con esta forma, podemos usar el metodo generico de avanzar sobre el rayo, en este caso solo vemos la parte de "arriba" del toro, que corresponde a la solucion positiva de la raiz.

GPURT: Funciones



Mientras que un rasterizador solo interpreta triangulos, el ray-tracing puede, en teoria, calcular la interseccion del rayo con lugares geometricos mas complicados. Por ejemplo una funcion de 2 variables z = f(x,y) que define una superfice en R3.
En lugar de armar una malla de triangulos a partir evaluar la funcion f en los vertices de control, se puede directamente resolver la ecuacion:

interseccion rayo - funcion
z = f(x,y)
x = rx + tDx
y = ry + tDy
z = rz + tDx

donde (rx,ry,rz) + t*(Dx,Dy,Dz) es la ecuacion del rayo.

Generalmente no hay forma analitica de resolver este sistema NO LINEAL, especialmente si la funcion es arbitraria y se comporta como una "caja negra". En lugar de intentar hayar una solucion analitica, se pueden usar metodos numericos.
La idea es que si z = f(x,y) y r(x,y) es la ecuacion del rayo, entonces se busca que
(1) f(x,y) - r(x,y) = 0
es decir el problema se reduce a encontrar la raiz de una funcion (NO LINEAL)
El metodo mas conocido para encontrar las raices de funciones arbitrarias es el metodo de Newton:

http://es.wikipedia.org/wiki/M%C3%A9todo_de_Newton

(ojo, que en nuestro caso se trata del Newton multivariable, y hay que usar la inversa de la matriz Jacobiana y no la derivada de la funcion!)

La eleccion del punto inicial en el metodo de Newton es un verdadero problema, y si no se elije un punto suficientemente cercano a la raiz, el metodo puede no converger.
Por este motivo, vamos a usar otra filosofia: partir del punto inicial del rayo, e incrementar un cierto dt, sobre la direccion, en cada paso evaluar (1) buscando un cambio de signo. Al encontrar el cambio de signo, aplicamos una busqueda binaria segun los cambios de signo hasta que la diferencia sea lo suficientemente pequeña.

Este metodo es bastante directo de implementar en la GPU. El mayor problema es determinar el dt para moverse sobre el rayo. Un dt demasiado grande puede hacer que el rayo pierda la interseccion y un dt demasiado chico haria que el algoritmo fuera demasiado lento. Usualmente se utiliza una malla de triangulos auxiliar, o un convex hull para hallar un punto cercano a la interseccion y partir de ahi. Eso no lo voy a hacer (porque es bastante tedioso de implementar, y ademas generar los triangulos va en contra de la filosofia de este ray-tracing).

void PS_RayTracing3( float4 Diffuse:COLOR0, float2 Tex : TEXCOORD0,
out float4 Color : COLOR0)
{
// Calculo la direccion del rayo u Obtengo la direccion del rayo
float3 D = normalize(g_vViewDir + g_vDx*(2*Tex.x-1) + g_vDy*(1-2*Tex.y)).rgb;
float near = 5;
float far = 40;
float dt = 0.25;
bool hit = false;
float3 pt = LF+D*near;
float dy = F(pt.x,pt.z)- pt.y;
int signo = sign(dy);
int cant = 0;
float t;
float color_actual = 0;
for(t=near;t<far&& cant<200 && !hit;t+=dt)
{
pt = LF + D*t;
dy = F(pt.x,pt.z)- pt.y;
if(sign(dy)!=signo)
hit = true;
else
++cant;
}

if(hit>0)
{
float t1 = t-dt;
float t0 = t1-dt;
for(int i=0;i<15;++i)
{
t = (t0+t1)/2;
pt = LF + D*t;
dy = F(pt.x,pt.z)- pt.y;
if(abs(dy)<0.001)
i = 10;
else
if(sign(dy)!=signo)
t1 = t;
else
t0 = t;
}

// aproximo la normal
float ds = 1;
float3 ipx = float3(ds,F(pt.x+ds,pt.z)-F(pt.x,pt.z),0);
float3 ipz = float3(0,F(pt.x,pt.z+ds)-F(pt.x,pt.z),ds);
float3 N = -normalize(cross(ipx,ipz));
float3 color_obj = tex2D( g_samDef, float2(pt.x,pt.z));
float3 g_vLightDir = normalize(g_vLightPos-pt);
float kd = saturate(dot(g_vLightDir,N)); // luz diffusa
float ks = saturate(dot(reflect(g_vLightDir,N), g_vViewDir));
ks = 0.75*pow(ks,5); // luz especular
color_actual = saturate((k_la + (1-k_la)*kd)*color_obj + ks);
}

Color.a = 1;
Color.rgb = color_actual;
}



Como se puede ver, hay algunos valores arbitrarios que hay que ir probando manualmente, como el dt=0.25 y la cantidad maxima de pasos una vez encontrado el cambio de signo (ahora en 15).
En la practica eso no representa mucho problema.

Por otra parte escribir la funcion f(x,y) es bastante sencillo, y no requiere introducir ningun otro dato adicional, (como la derivada).

float F(float x,float y)
{
return 35*exp(-2*(0.008*(x*x+y*y)))*sin(0.008*(x*x-y*y));
}







El algortimo es estable siempre y cuando la funcion sea lo suficiemente suave. Mientras que teoricamente la continuidad y Bolzano aseguran la existencia de una solucion, el dt arbitrario y la estabilidad numerica del algoritmo requieren mas que una continuidad C1. Si la funcion cambia demasiado rapido en un pequeño entorno del rayo, este puede perder la interseccion: como en este caso, con la funcion 1/(x+y) (que dicho sea de paso no es continua!)




La combinacion de discontinuidades tanto en la funcion como en su derivada, sumado a los errores de precision propios del algoritmo generan inestabilidades que recuerdan a los atractores fractales. En el fondo el concepto es bastante similar: Imaginen la funcion sin(1/x) para todo x!=0. Aunque la funcion es continua (salvo para x=0), en un entorno de x=0, cambia de signo cada vez mas veces. Por mas precision (dt cada vez mas pequeño) es imposible que el mismo dt sirva para todos los puntos en dicho entorno. La cantidad de cambios de signo tiende a infinito: el algoritmo se hace inestable.
En resumen, hay puntos p0 en los que el algoritmo converge pero en un punto infinitamente cercano a p0, sea p0 + epsilon, el algoritmo diverge, generando imagenes caoticas, como estas :



viernes, 19 de noviembre de 2010

GPU-RT: Reflexiones

Una de las grandes ventajas del ray-tracing con respecto al rasterizado es la posibilidad de calcular reflexiones. La reflexion es simplemente calcular un rayo mas: hay que hacer rebotar el rayo, y seguir haciendo el ray-tracing de forma recursiva tantas veces como el nivel de detalle que se requiera.

Sin embargo, en la GPU eso plantea unos cuantos problemas, ya que no se podemos llamar al PS dentro del PS, y tenemos que ir generando la imagen en etapas sucesivas.

Para resolverlo, vamos a plantear un sistema de "retroalimentacion". Primero, como siempre se dibujan los triangulos a efectos de llamar al PS para cada pixel. Pero la salida se redirige a 3 texturas simultaneas (se requiere que la placa soporte multiples renders target, si no lo tiene igual tambien se puede hacer, pero se requieren todavia mas pasos).
El primer RenderTarget sera el color buffer pp dicho. El segundo el rayo reflejado, y el tercero el punto de interseccion.
En el segundo paso, estas 3 texturas se usan como datos de entrada. El rayo reflejado se transforma en la direccion del rayo pp dicha. Y el colorbuffer anterior sirve para hacer el blending. El punto de interseccion se usa para calculos internos.

void PS_RayTracing( float4 Diffuse:COLOR0, float2 Tex : TEXCOORD0,
out float4 Color : COLOR0 ,out float4 RayDir : COLOR1,out float4 oIp : COLOR2 )
{
// Calculo la direccion del rayo u Obtengo la direccion del rayo
float3 D = nro_R==0?
normalize(g_vViewDir + g_vDx*(2*Tex.x-1) + g_vDy*(1-2*Tex.y)).rgb:
tex2D( g_samRayDir, Tex);

......
......


Como se ve en el codigo, la primera vez (nro de recursion = nro_R = 0) se calcula el rayo, pero las siguientes se obtiene el rayo de la salida anterior.
Entonces, cuando se requiere que el algoritmo siga calculando rayos, por ejemplo el caso de un espejo, hay que generar el rayo reflejado en salida:

if(objsel!=-1)
// Calculo la direccion del rayo reflejado
D = reflect(D,N);
else
// el 2 indica que no hay mas rayos que trazar
D.x = 2;
.....
.....
// out RayDir
RayDir.xyz = D;
RayDir.w = 0;




Al final del PS se genera la direccion del rayo reflejado en salida, que sera guardado en una textura auxiliar, para usar como dato de entrada en el siguiente paso.








piso reflejando las esferas, que a su vez reflejan el piso y si mismas. Hay 2 niveles de reflexion.



misma escena con 3 niveles de recursion:


misma escena con 6 niveles de recursion:


Ahora, agregando colores a las esferas + luz especular, segun estas formulas:

float kd = saturate(dot(g_vLightDir,N));        // luz diffusa
float ks = saturate(dot(reflect(g_vLightDir,N), normalize(Ip-LF)));
ks = 0.7*pow(ks,5); // luz esspecular
color_actual = saturate((k_la + (1-k_la)*kd)*color_obj + ks);






En cada paso, los resultados guardados en las texturas del paso anterior, se combinan con los actuales con una operacion del blend. Para el blend se puede usar un factor de alfa constante, o bien, que ese valor se vaya haciendo cada vez mas pequeño en cada reflexion, motivado por el hecho que los rayos reflejados pierden energia en cada reflexion, por eso las imagenes se ven mas tenues en las sucesivas recursiones.


if(nro_R==0)
Color.rgb = color_actual;
else
{
//float alfa = pow(0.5,nro_R);
float alfa = 0.5;
Color.rgb = alfa*color_actual + (1-alfa)*tex2D( g_samColorBuf, Tex);
}

Ejemplo con alfa = 0.5


GPU-RT : Mipmapping

La imagen anterior tenia un problema, en los bordes de la esfera se ve un color incorrecto, que parece un problema de aliasing:




Sin embargo el defecto no tiene nada que ver con el aliasing. De hecho, si se reemplaza el color de las esferas por un color fijo y el piso tambien, se puede observar que no hay ningun problema en el contorno, y la precision es la del pixel.




El problema es mucho mas complicado y esta relacionado con el mipmapping de las texturas. Cuando se carga una textura el directX puede generar una serie de texturas reduciendo el tamaño por la mitad en el ancho y en el alto, y generando una cadena de texturas, que facilitan el "texture filtering". Asi para una textura de 64x64, se generan otras de 32x32, 16x16...hasta una de un pixel.

El pipeline automaticamente selecciona la textura correspondiente haciendo uso de lo que se llama la derivada parcial de la coordenada.

Ejemplo: del manual del direct:
ddx - HLSL
Returns the partial derivative of x with respect to the screen-space x coordinate.

The x is supposed to be an interpolated parameter, for example a texture coordinate. The ddx function indeed returns the rate of change of it's parameter across the current pixel and the adjacent one. This works, because hardware that supports the derivative operations will always rasterize at least a 2x2 grid of pixels at the same time.

El problema que tenemos es que en verdad estamos dibujando 2 triangulos! Y cuando se llega al borde con otro objeto la grilla interna del hardware va a tener pixels que corresponden a coordenadas de textura completamente diferentes (inclusive de diferentes texturas). La ddx le da cualquier cosa, y usa un mimmap incorrecto, generando ese problema en el borde. Usualmente le da un valor enorme, y utiliza el ultimo mimap, que corresponde al promedido de todos los puntos de la textura. (Este hecho es el que me ayudo a darme cuenta lo que pasaba), porque en el borde me ponia el promedio de la textura!

Desafortundamente, resolverlo no es tan facil. Para hacerlo bien habria que calcular 2 rayos adicionales, moviendose un pixel en x y en y respectivamente, de tal forma de determinar a que punto de la textura van a parar, y ahi determinar cuantos texels se requieren promediar, para buscar el bitmap mas adecuado.

Esto quiza sugiere modificar el algoritmo para que trabaje con 4 rayos al mismo tiempo, pero por el momento, una forma facil de resolverlo es usar la funcion
tex2Dlod en lugar de tex2. Esta funcion recibe como parametro tambien el LOD de la textura. El mismo se puede aproximar usando la distancia desde la camara al punto de interseccion, basicamente cuanto mas lejos este de la camara menos nivel de detalle utiliza, o sea mas grande el parametro que se le pasa a tex2Dlod.
Otra vez, no hay una forma de calcularlo rapido, y hay que tocarlo a "dedo".

Ejemplo:

// interpolo la textura en el triangulo                
float2 tex = txA*(1-bc_b-bc_g) + txB*bc_b + txC*bc_g;
//float3 color_obj = tex2D( g_samDef, tex);
float3 color_obj = tex2Dlod( g_samDef, float4(tex.x,tex.y,0,R/35));


R es la distancia desde la camara hasta el pto de inteseccion, pero el 35 surge de la distancia maxima visible y la el hecho que la textura tenia 8 mipmapings, o sea esta a "ojo". En general una funcion lineal como R*K resuelve el problema en casi toda la imagen, menos en lo que esta muy lejos, como se puede ver en la siguiente imagen, que los contornos de la esfera salen bien, pero la parte de atras del piso, no tiene mucha resolucion.

GPU-RT: Interseccion con el Triangulo

La interseccion del rayo con el triangulo (paradojicamente) es mas complicada que con la esfera, y ademas trae aparejado un nuevo problema: el texture mapping.

Probablemente este es uno de los motivos por el cuales el rasterizer es mas eficiente que el ray-tracing. El triangulo desde el pto de vista del rasterizador, solo tiene que calcular 3 vertices y despues interpola todo. Pero para el raytracing son todos rayos independientes en los cuales tiene que calcular intersecciones.

El mayor problema con la interseccion entre el rayo el triangulo es que algebraicamente podemos hayar facilmente la interseccion entre la recta y el plano. Pero luego tenemos que verificar que el punto este dentro del triangulo.

Afortunadamente, este tema esta sumamente estudiado, y usualmente se usan las coordeandas baricentricas para expresar la posicion de un punto dentro de un triangulo.

(solo la parte que habla del triangulo)
http://en.wikipedia.org/wiki/Barycentric_coordinates_(mathematics)

Pero se los resumo, las coordenas baricentricas de un punto son P, con respecto a un triangulo de vertices A,B y C son, b y g tal que :

P = A*(1-b-g) + B*b + C*g;
con b>=0, g>=0

Y tienen la propiedad que si el punto esta adentro del triangulo se cumple que b+g<=1, de hecho si b+g = 1 esta justo en la frontera del triangulo.

Las coordenadas baricentricas nos resuelven 2 problemas al mismo tiempo: primero si el punto esta adentro del triangulo, y ademas despues nos sirven para interpolar cualquier valor definido en los 3 vertices, para cualquier otro punto del triangulo. Asi es como funciona intermanente la GPU. Lamentablemente nosotros no estamos usando el pipeline, con lo cual lo tenemos que reprogramar a manopla:

Haciendo un poco de cuentas, para calcular la interseccion y las baricentricas necesitamos resolver unos 3 determinantes, afortunadamente la funcion de determinante esta acelerada por hard, y hay una funcion del hlsl.

El codigo queda asi:

// Triangulo
//------------------------------------------
float3 A = tex2D( g_samObj, float2(j+0.5,i+0.5)/OBJ_SIZE);
float3 B = tex2D( g_samObj2, float2(j+0.5,i+0.5)/OBJ_SIZE);
float3 C = tex2D( g_samObj3, float2(j+0.5,i+0.5)/OBJ_SIZE);
float3x3 M = float3x3( B-A,C-A,-D);
float det = determinant(M);
if(det!=0)
{
float3x3 M3 = float3x3( B-A,C-A,LF-A);
float b,g,t;
t = determinant(M3)/det;
if(t>0 && t<R)
{
float3x3 M1 = float3x3( LF-A,C-A,-D);
b = determinant(M1)/det;
if(b>=0)
{
float3x3 M2 = float3x3( B-A,LF-A,-D);
g = determinant(M2)/det;
if(g>=0 && b+g<=1)
{
// nota: Ip = A*(1-b-g) + B*b + C*g;
// coordenadas baricentricas: a , b y g
// (a+b+g = 1)
R = t;
Ip = LF+D*t; // punto de interseccion
objsel = t;
i_sel = i;
j_sel = j;
tipo_sel = 1;
bc_b = b;
bc_g = g;
// calculo la normal
N = -normalize(cross(B-A,C-A));
}
}
}
}




Mas adelante (en el "pipeline" del raytracing), cuando se necesita calcular el color del pixel, lo unico que se tienen son las coordenadas baricentricas donde el rayo intesecto al triangulo.
Para relacionarlo con la textura, necesitamos mas informacion, asi que seguimos agregando datos a los objetos. El caso del triangulo es uno de los que mas datos requiere: 3 posiciones x 3 vertices + 2 coordenadas uv x 3 vertices, de momento tenemos 4 texturas formato D3DFMT_A32B32G32R32F (llamemosla g_pTexObj,g_pTexObj2, etc) + una textura D3DFMT_X8R8G8B8, (llamada g_pTexObjT)
La textura g_pTexObjT de momento solo usa el canal R para el tipo de objeto (0==esfera, 1==triangulo, etc)
Las texturas g_pTexObj guardan el resto de la info (para el caso del triangulo)
g_pTexObj.x -> coordenada x del primer vertice
g_pTexObj.y -> coordenada y del primer vertice
g_pTexObj.z -> z del primer vertice
g_pTexObj.w -> coordenada u de textura primer vertice

g_pTexObj2.x -> coordenada x del segundo vertice vertice
etc etc..

Las coordendas v de la textura (que no entraban) hubo que meterlas en
g_pTexObj.x = coordenada de textura V del primer vertice
g_pTexObj.y, g_pTexObj.z idem segundo y tercer vertice.

Bueno, cuando llega el momento de tomar el color del pixel de la textura, se acceden a todas estas texturas, y si interpola, usando las baricentricas:

// nota: Ip = A*(1-b-g) + B*b + C*g;
// intersecto en el triangulo
// tomo los datos del objeto
float4 A = tex2D( g_samObj, float2(j_sel+0.5,i_sel+0.5)/OBJ_SIZE);
float4 B = tex2D( g_samObj2, float2(j_sel+0.5,i_sel+0.5)/OBJ_SIZE);
float4 C = tex2D( g_samObj3, float2(j_sel+0.5,i_sel+0.5)/OBJ_SIZE);
float4 D = tex2D( g_samObj4, float2(j_sel+0.5,i_sel+0.5)/OBJ_SIZE);
// tomo las coordenadas uv de cada punto
float2 txA = float2(A.w,D.x);
float2 txB = float2(B.w,D.y);
float2 txC = float2(C.w,D.z);
// interpolo la textura en el triangulo
float2 tex = txA*(1-bc_b-bc_g) + txB*bc_b + txC*bc_g;
float3 color_obj = tex2D( g_samDef, tex);
float kd = dot(N,g_vLightDir); // luz diffusa
color_actual = (k_la + (1-k_la)*kd)*color_obj;




Bueno...vamos a tratar de ir llegando a la imagen del wikipedia....sacando las diferencias!!!

GPU-RayTracing

Alla por mediados de los 90, hice mi primer intento de escribir un ray-tracing. Estaba tratando de obtener imagenes “realistas” de muebles, para un software de cálculo. Algunos años antes de tener acceso a internet, habia leido sobre el algoritmo en un libro de trucos graficos 3d en GWBasic.
El raytracing es uno de esos algoritmos que en la teoria son tan fáciles que dan miedo:

http://es.wikipedia.org/wiki/Raytracing

y cuando uno los quiere implementar en la practica....si no se trabaja mucho sobre el tema, o la calidad es pobre o tarda 2 horas en dibujar una imagen. (o las 2 cosas)

El tema me hace acordar a esos algortimos para calcular numeros primos, que pueden ir desde unos cuantos for anidados y que tarde lo que tarde, hasta series de Rieman que no las entiende ni Riemann.


El raytracer que escribi se basaba en soportar la mayor cantidad de primitivas: esferas, elipses, conos, toros, y otras superficies cuarticas. En ese momento los triangulos no me parecian figuras muy atractivas, y le dedique mas esfuerzo a solucionar las ecuaciones de intereccion entre el rayo y superficies mas complicadas como toros.

LOAD TEXTURE piso4
LOAD TEXTURE piso3
LOAD TEXTURE marmol5
LOAD TEXTURE hela
WORLD <27000,27000,5000>
SHADOWS OFF
LOOK FROM = <3000,4200,2100>
LOOK AT = <0,0,500>
VUP = <0,0,1>
FIELD OF VIEW = 0.9
LIGHT = <3000,3200,2000>
SCREEN = 600,400
DEFINE OBJECT COPA <100,100,100>
SURFACE = 0.5,0.9,0
COLOR = 128,128,192
TEXTURE = -1,0,0
OBJECT CIRCLE <-50,-50,0> , <50,-50,0>, <-50,50,0>
OBJECT CILINDER <0,0,0> , <0,0,50> , 10
OBJECT ELIPSOID <0,0,100> , <50,50,75> , <-75,0>
END OBJECT

............
............

OBJECT SILLA <1200,1500,0> , <300,300,1000> , <0,0,0>
OBJECT APARADOR <1600,0,1300>, <300,1000,700>, <90,0,0>
OBJECT APARADOR <0,1500,1300>, <300,1000,700>
OBJECT APARADOR <0,0,1300>, <300,900,700>
OBJECT EXTRACTOR <1600,10,1300> , <500,400,800>
OBJECT MICROONDAS <300,10,800>, <500,400,400> , <0,0,0>
END OBJECT

OBJECT COCINA <0,0,0>, <5000,5000,2300>


Ray-tracing viejo:






Desafortunadamente en el mundo real (especialemente el de los muebles) no hay muchos toros, y si hay triangulos (o rectangulos) por todos lados.

Unos cuantos años mas tarde quiero retomar el tema pero desde otro punto de vista. La idea es implementar el ray tracing en la GPU, en lo que seria un ejemplo de usar la GPU para proposito general. (Vean CUDA en la pagina de nvidia)

Si bien todavia el raytracing esta lejos de competir con un rasterizer en aplicaciones de real time, hay algunas ventajas interesantes, especialmente con algunas geometrías.
Por ejemplo el caso de una esfera: para representarla en un raytracing se puede usar sola la posicion y el radio. Y luego, la precision de la misma es “ilimitada”, en el sentido que por mas que nos acerquemos la precision siempre va a ser la misma.
Para lograr lo mismo con un rasterizador, habria que transformar la esfera en triangulos, cientos de ellos. Y por mas triangulos que usemos, siempre que nos acerquemos lo suficiente se veran los triangulos.


Paradojicamente lo que queremos hacer es implementar el ray-tracing en la GPU. Que ventaja pueda tener? Mientras que no se si va a ser mejor o peor que el mismo algoritmo en la CPU hasta que no lo pruebe efectivamente (y dudo que lo haga) si hay 2 cuestiones por las cuales la GPU puede llegar a ser convieniente:
- Por el alto grado de paralismo que se puede lograr, siendo que cada rayo es completamente independiente de todos los demas, es un proceso ideal para correr como un shader en la gpu que los va a ejecutar varios al mismo tiempo.
- Una vez que se tiene el pto de intereccion se puede usar la funcion del shader tex2d para texturizar y evitarse programar los filtros de texturas por software.



Estrategia.

La idea es meter el algoritmo de interseccion entre el rayo y todos los objetos de la escena en un Pixel Shader. El PS recibira como parametro la direccion del rayo, y tendra que devolver el objeto mas cercano con el cual intersecta.

Como almacenamos los datos de los objetos.
El PS tendra que acceder a todos los objetos para ir calculando intersecciones, con lo cual vamos a usar una (o muchas) texturas para guardar la info de los mismos. Para la esfera necistamos 3 float para la pos y un float para el radio. Pero tambien necesitamos el tipo de objeto, y algunos objetos tienen mas datos, paradojicamente los triangulos tienen 3 puntos, lo que da 9 float.

Ejemplo,
creo una textura de un cierto tamaño suficientemente grande, de formato float.

g_pd3dDevice->CreateTexture( OBJ_SIZE,OBJ_SIZE,
1, D3DUSAGE_DYNAMIC,
D3DFMT_A32B32G32R32F,
D3DPOOL_DEFAULT,
&g_pTexObj,NULL);
g_pEffect->SetTexture( "g_txObj", g_pTexObj);




Luego, tengo que inicializarla con los datos de los objetos pp dichos:

D3DLOCKED_RECT lr;
g_pTexObj->LockRect( 0, &lr, NULL, 0);
BYTE *bytes = (BYTE *)lr.pBits;
FLOAT *texel = (FLOAT *)(bytes+t);
// formato del texel (4 bytes por canal)
// | R | G | B | A | = 128 bits
texel[0] = 3; // pos X
texel[1] = 2; // Pos y
texel[2] = 0; // Pos Z
texel[3] = 0.5; // Radio
// etc etc etc




Como hacemos que se ejecute el PS.
Como hay que trazar un rayo desde el punto de vista hacia cada pixel de pantalla, eso nos sugiere dibujar 2 triangulos que ocupen toda la pantalla, y usar la pos. en pantalla para calcular la direccion del rayo, y asi se llamara el PS para cada uno de ellos.


El PS (resumido) queda asi:

void PS_RayTracing( float4 Diffuse:COLOR0, float2 Tex : TEXCOORD0,
out float4 Color : COLOR0)
{
// Calculo la direccion del rayo u Obtengo la direccion del rayo
float3 D = normalize(g_vViewDir + g_vDx*(2*Tex.x-1) + g_vDy*(1-*Tex.y)).rgb;
float R = 100000;
float3 Ip = 0;
int objsel = -1;

// recorro todos los objetos
for(int t=0;t<cant_obj;++t)
{
// tomo los datos de la esfera
float4 obj = tex2D( g_samObj, loat2(j+0.5,i+0.5)/OBJ_SIZE);
float3 esfera = obj.xyz;
float radio = obj.w;

// verifico la interseccion entre el rayo y la esfera
// si intersecta, y en una distancia menor a la actual (R)
// calculo el pto de interseccion
Ip = LF+D*t1; // punto de interseccion
R = tl;
objsel = t;

// paso al siguiente objeto
}

// si intersecto con algun objeto: tomo el color correspondiente
// al punto donde intersecto:
if(objsel!=-1)
// intersecto en la esfera objsel
// calculo mapping u-v de la esfera etc etc
// tomo el color de la textura del objeto
color_obj = tex2D( g_samDef, float2(u,v));
else
color_obj = color_fondo; // no toca contra nada


// devuelve el color del pixel
Color.rgb = color_obj;
Color.a = 1;
}







Como se puede ver, el algoritmo es bastante sencillo. Sin embargo la logica esta “forzada” para adaptarse a la GPU (como siempre que se quiere programar algo en la GPU que no sea dibujar triangulos).

Algunos de los problemas:
Si bien eso puede variar de GPU en GPU, (lo importante es el concepto), haciendo pruebas hay un cierto limite en la cantidad de objetos que se pueden evaluar en el PS. Por ejemplo con unos 100 objetos no hay problema, pero si fuesen 1000 objetos no soporta esa cantidad de lineas el PS (Los for....se explotan en varias lineas, no es como la CPU)
Eso nos sugiere que, al igual que en la CPU, hay que implementar una estructura jerarquica para organizar los objetos en el espacio, de tal forma de mantener la cantidad de intersecciones que se evaluan en cada PS en no mas de 100.
A diferencia de la CPU, aca tenemos que ver meter esta estructura en texturas, y dividir el proceso en pasos intermedios e ir almacenando los resultados parciales en texturas auxiliares.
Pero eso, para la proxima....


Poniendo todo junto:
En la primera prueba, vamos a usar solo 100 objetos, todos esferas, a modo de prueba.