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.





domingo, 31 de octubre de 2010

Shadow-Texture Map

Este post esta relacionado con el de Shadows Map anterior. La idea surgió haciendo modificaciones al shadow map para que en lugar de proyectar una sombra, proyecte una luz o una textura. (ojo que no es un light map), por ejemplo una persiana que proyecta la luz (y la sombra) sobre el piso.
Esto mismo se podría hacer generando la geometría de la persiana y ubicando una luz atras de la misma, y usando el shadow map standard, pero eso es muy costoso, en el sentido que la geometría suele ser bastante complicada. (En verdad ni lo probé).
Lo que si probé, es modificar la textura del shadow map a mano, simulando que existe una geometría, y eso si funciona relativamente rápido. Por ejemplo, si escribo un texto en el shadow map (una vez que lo tengo en memoria ya calculado, usando el shadow map en 2d), el efecto es que se proyecta el texto como sombra en la escena. Sin embargo después de analizarlo más a fondo, encontré que ni hace falta generar el shadow map, y la cosa es mucho más sencilla.
Vamos de cero:
En el SM standard, el problema central radica en obtener las coordenadas del shadowmap a partir de las coordenadas del pixel actual. En otras palabras, lo que quiero saber es que coordenadas tendría el pixel actual desde el punto de vista de luz, y no de la camara.

Si partimos de la posición iPos de un vértice en object space y en lugar de aplicar la transformación de la camara, aplicamos la de la luz:

vPosLight = mul( iPos, g_mLightWorldViewProj );

g_mLightWorldViewProj es lo mismo que el g_mWorldViewProj pero haciendo de cuenta que la camara esta posicionada en la luz, y mirando hacia la direccion arbitraria a la que apunta el spot de la luz pp dicha.


Ahora, en el shadowmap standard, habria que usar las coordenadas vPosLight, para buscar en el shadowmap (previamente calculado) y verificar si el punto esta iluminado o no. Sin embargo lo que vamos a hacer es mucho mas simple:
1- No hay ningun shadow map si no una textura previamente cargada, que es la que queremos proyectar
2- En el VertexShader generamos vPosLigth como siempre: vPosLight = mul( iPos, g_mLightWorldViewProj );
3- En el PixelShader, usamos vPosLight para calcular la pos. en la textura a proyectar, y si esta dentro de los limites de la misma (x,y entre 0 y 1), tomamos el color de la textura.

Ejemplo:

float4 PixScene2( float2 Tex : TEXCOORD0,
float4 vPos : TEXCOORD1,
float3 vNormal : TEXCOORD2,
float4 vPosLight : TEXCOORD3 ) : COLOR
{
float4 color1 = 0;
float4 color2 = tex2D( g_samScene, Tex );
if(vPosLight.w>=0)
{
float2 pos = 0.5 * vPosLight.xy / vPosLight.w + float2( 0.5, 0.5 );
pos.y = 1.0f - pos.y;
if(pos.x>=0 && pos.x<=1 && pos.y>=0 && pos.y<=1)
color1 = tex2D( g_samProjTex, pos);
}

return lerp(color1,color2,0.5);
}










Esta linea sirve para limitar la proyeccion a la parte que esta delante del nearplane
if(pos.x>=0 && pos.x<=1 && pos.y>=0 && pos.y<=1)
De lo contrario saldria asi:





Por otra parte, la matematica involucrada en la proyeccion genera 2 imagenes, una inversa de la otra, cuando las coordenadas homogeneas son ambas negativas, es decir
tanto z, como w son negativas, luego Zp = z/w >0, pero estamos del otro lado de la "camara", ejemplo:




Para eliminar ese artificio, simplemente se puede poner:
if(vPosLight.w>=0)

Por ultimo, habría que combinarlo con el ShadowMap estandard, para evitar proyectar la textura en los lugares donde no podria verse, porque esta oculta por otro objeto.

lunes, 6 de septiembre de 2010

El Plano Proyectivo

Consideremos un punto P=(1,0) en el plano, la recta real en y=0, y un círculo que pasa por el origen y por P.

Graficamente:



Para cada punto de la recta real , se puede trazar una recta hasta el punto P, que intersecta al círculo en un único punto. Esta es la demostración standard, que el intervalo abierto (-1,1) tiene la misma cantidad de puntos que la recta real. Primero se demuesta que el intervalo (-1,1) se deforma hasta formar una circunferencia a la que le falta un punto (por eso el intervalo es abierto y no cerrado). Luego, la recta que une los puntos de la circunferencia con la recta real, define una función biyectiva que es lo que se necesita para probar que ámbos conjuntos tienen la misma cantidad de elementos.

La función se puede escribir como
(1) x1 = x1p / h
para todo x1p que pertenece a la circunferencia, salvo el punto p = (1,0), que es donde h=0 y el cociente no existe.

El hecho que la circunferencia tenga un punto menos que la recta real, a los efectos de la cardinalidad no tiene importancia, y se demuestra facilmente que un número finito o infinito numerable de elementos no modifica la cardinalidad de un conjunto infinito.
Sin embargo si genera una diferencia importante en el aspecto topológico, la recta real es equivalente a la circunferencia a la que le falta un punto, no la circunferencia completa.

Interpretación.
Si queremos representar en la recta real todos los ángulos posibles, o sea todas las rectas que pasan por P, nos estaría faltando una recta, que es la paralela a la recta real. Si se agrega el punto que falta en la circunferencia y se lo relaciona con dicha recta, (intuitivamente seríaa agregar el punto del infinito), podemos establecer que el conjunto de la recta real unión el punto del infinito es topologicamente una circunferencia completa.
A este conjunto se lo llama la recta proyectiva, porque esta compuesto de todas las rectas que pasan por un punto, o dicho de otra manera, por la recta real más el punto del infinito.
Los 2 puntos en el infinito se toman como el mismo punto en realidad, ya que corresponden a una sola recta, la dimensión de los puntos del infinito es cero, mientras que la dimensión de la recta es 1.


Plano proyectivo.
Pasando a 2 dimensiones, y con el mismo esquema de demostración, se tiene que el plano es equivalente a una esfera (superficie, usualmente se llama S2) a la que le falta un punto, o tambien, a una semiesfera abierta (que no incluye el borde). La frontera (el ecuador) esta relacionada con los puntos en el infinito. Se puede definir al plano proyectivo como el conjunto de todas los puntos del plano más los puntos en el infinito. En este caso los puntos del infinito son equivalentes a toda una circunferencia, y su dimensión es 1 , mientras que el plano tiene dimensión 2. De esta forma se puede definir recursivamente para un hiperplano:
El espacio proyectivo, se define de la misma manera, y los puntos del infinito son un plano y asi sucesivamente.


Coordenadas Homogeneas.
Para describir un punto en la recta proyectiva si bien se podría usar un sólo elemento ya que su dimensión es 1, (el agregarle un punto no cambia la dimensión del espacio), usualmente se usan 2 elementos, relacionados por la formula (1), que la escribimos de nuevo:

(1) x1 = x1p / h

Q = (x1p,h) cada elemento del plano proyectivo se representa por 2 coordenadas, pero que no son independientes si no que estan relacionadas por la fórmula (1).

Esto llevado al plano proyectivo, y usando la noticion habitual, se tiene
Q = (xp,yp,w)
donde x = xp/w y y=yp/w y (w!=0)

Esta forma de notación se llama coordeandas homogeneas o w-homogeneas. Se interpreta que el punto esta siendo proyectado en el plano w=w0, usualmente para puntos en el plano se utiliza w=1, con lo cual se tiene:
Q = (xp,yp,1)

Y tiene la ventaja que es posible escribir los puntos en el infinito, que estan asociados a w=0,
Q = (xp,yp,0)

Siempre se hace énfasis en que xp/0 no existe, y que no se tiene que tomar como una división que da infinito. Uno intuitivamente interpreta que xp/0 es un símbolo que representa un punto en el infinito y no una operación de división.

Algunas consecuencias de esta notación son
- cada punto tiene infinitas representaciones,
- que se puede escribir los puntos del infinito.
- 2 rectas siempre se cortan en un punto.

Estas coordeandas fueron introducidas por Möebius hace como 200 años pero se usan hoy en dia en la mayor parte de las librerias gráficas.

Que esta estudiando Möebius y que tiene que ver con el plano proyectivo?

Imaginen que queremos representar todas las rectas que pasan por un punto P en el espacio, (ese conjunto es el plano proyectivo claro). La forma mas facil es ubicar el punto P en el centro de una esfera (superficie), y las rectas que pasan por el centro, intersectadas con la esfera generan 2 puntos: p1,p2
Sin embaro, no hace falta usar 2 puntos, si no que con uno solo ya alcanza, porque la recta siempre pasa por el origen, con lo cual con una semiesfera es suficiente:
Todas las rectas están representandas por un punto en el interior de la semiesfera S, salvo (las del infinito) que van a estar representadas por 2 puntos en la frontera de S, o sea en el ecuador. Si arbitrariamente tomamos uno de los 2 puntos, (no hacen falta los 2), y si tomamos una recta que se va moviendo, siempre representada por un punto en S, cuando llega al ecuador (la frontera de S), súbitamente desaparece y re-aparece por el otro lado en la antípoda. (Eso en el momento que se hace invisible el primer extremo de la recta, a la vez que se hace visible el segundo, por el otro lado)



Como podemos “reorganizar” la superficie para que el movimiento de este punto virtual sea continuo?

Que pasa si intentamos “cocer “ la semiesfera S por el borde, uniendo cada par de antipodas? Bueno, eso no se puede hacer, intenten con una tela, o se puede demostrar con topología, (si se pudiese la esfera sin un punto seria equivalente a la que esta completa y eso no es cierto).

Möebius (o vaya a saber quien, la verdad que no lo se, pero me imagino que habra sido él), propuso lo siguiente:

Tomamos un punto p1 en el ecuador, y recortamos un entorno sobre el, y sobre su opuesto p2. Esas 2 superfices las unimos invertidas por el borde que antes era el ecuador , y con eso obtenemos un disco, el pto p1 = p2 = centro del disco, y con esa operación conectamos los puntos p1 y p2, y todo el entorno de ellos en el ecuador.
Lo quedo es esto:



Tomemos otro punto q1 en lo que quedo del ecuador, su opuesto q2, esta del otro lado, pero “opuesto”. Quiero conectar esos 2 puntos, (ya que son el mismo elemento), para ello tengo que doblar de un extremo y unir invirtiendo, lo que queda es una banda de Möebius, (que casualidad no? )




Ahora, la banda es una superficie no orientable: tiene una sola cara y un solo borde cerrado. Si el borde es una curva cerreda, es equivalente a una circunferencia. Si intentamos doblar y plegar la banda para que nos quede ese círculo no vamos a poder en 3 dimensiones, se necesitan 4 dimensiones para que la banda no se intersecte a si misma, pero eso no quiere decir que no se pueda imaginar. Sigan el borde de una banda de Möebius y van a verificar que es una curva cerrada, despues olvidense de la superficie, imaginense solo una curva cerrada y vean que se puede deformar hasta convertirse en una circunferencia.

Bueno, una vez que se tiene el círculo, ese borde corresponde a la frontera de los 2 entornos que habiamos retirado originalmente, y que unimos para formar un disco. Si pegamos ese disco en el borde, quedan conectados todos los puntos.
Con un poco de esfuerzo se puede demostrar que esta construcción respeta la "estructura" original (es un morfismo especial), si por ejemplo estudiamos el movimiento de una recta en la esfera original, representandola como el movimiento del punto en el plano proyectivo, el movimiento es continuo en ambos casos.

La banda de Möebius unida por el borde con un disco, se llama en topologia el plano proyectivo.

domingo, 5 de septiembre de 2010

Refrito de Matemáticas - Parte A

Esto esta reflotado de un TP de cuando estaba cursando análisis funcional y álgebra lineal en la UBA. Son todas cuestiones de matemáticas que tienen que ver con los gráficos por computadora.

El producto interno.
Motivación Geométrica.

Consideremos 2 vectores en el espacio euclídeo A y B, y el producto interno definido
A . B = A B cos O
Donde O = ángulo que forman, graficamente:

A cos O representa la proyeccion escalar de A sobre el vector B

Se puede demostrar que en R3 es equivalente al producto matricial:

Bx
A.B = [Ax Ay Az] By = AxBx + AyBy + AzBz
Bz


A veces se suele notar el producto interno entre 2 vectores v y w, como , especialmente cuando se trata de otros espacios que no sean Rn

Si consideramos dos vectores v y w en R3, w =1 (unitario) y el subespacio W (de dimensión en 1) generado por w, es decir W = {k*w}
Se puede demostrar que el vector
p = <v,w>*w
que es la proyección de v sobre w, es la mejor representación de v en W. Es decir no existe otro elemento en W que aproxime mejor a v.

Al estudiar el producto interno en el espacio euclídeo (recta, plano y espacio), los matemáticos indentificaron cuales eran las mínimas propiedades necesarias para demostrar todos los teoremas y asi generalizaron el concepto del producto interno a espacios mas abstractos, como los espacios de funciones. Finalmente se abstraen de la fórmula del producto interno, y formalmente, se define como una aplicación de VxV en K: (V es el espacio Vectorial y K es el cuerpo sobre el que esta definido. )


<.,.> VxV -> K

Ahora <.,.> tiene que cumplir:

<ax + by,z> = a<x,z> + b<y,z> (linearidad)

<x,y> = conj <y,x> (hermeticidad)

<x,y>=0 y <x,y>=0 <=> x=0 (definida positiva)





Es decir cualquier aplicación que cumpla esas condiciones es un producto interno, y son válidas todas las afirmaciones que se demostraron para Rn
Es facil demostrar que A.B = AxBx + AyBy + AzBz es un producto interno en R3, lo mismo que A.B = A1B1 + A2B2 + ....+AnBn tambien es un producto interno en Rn
(es el primer ejercicio que se ve en álgebra luego de definir el producto interno)

Usualmente esta definición abstracta de producto interno es el punto de partida del cual se deducen todas las demas propiedades. Por ejemplo, el angulo entre 2 vectores, se define a partir del producto interno como

cos O = <v,w>/(v w )

Esto permite definir un ángulo para vectores en especios vectoriales no tan convencionales como Rn. La verdad que la gracia de tanta abstracción se aprecia cuando se trabajan con otro tipo de espacios, como los espacios de funciones, donde un vector es una función, y gracias a estas abstracciones se puede definir el “ángulo” entre 2 funciones.

Bases.
De algebra lineal sabemos que todo vector se puede expresar como una combinación lineal de elementos de la base. Esa combinación es única, y sus coeficientes se llaman coordenadas.
Ahora, si la base es ortonormal, o sea si los vectores son ortogonales entre si, y unitarios, se cumple que :


B = { b1,b2...bn} Base ortonormal

v = <v,b1>b1 + <v,b2>b2 +...+ <v,bn>bn



Demostración. (resumida, para ver la idea)


por ser B base

v = a1b1 + a2b2 + ...+ anbn

lo que queremos probar es que

vk = <v,bk>para todo k

<v,bk>= <sum>x las prop de linealidad de <> para la suma y para el producto:

= sum ai <bi,bk>= ak porque <bi,bk>= 0 si i==k, o =1 si i=k

como queriamos probar,



Quería demostrar este teorema para hacer hincapie en que sólo se usa la definición del producto interno y la ortonormalidad de la base, y no la fórmula concreta para calcularlo.
Esta ecuación a la conoce tambien como expansión de Fourier, por algo que despues vemos.
Ahora, si consideramos la matriz T = [b1,b2,...bn], es decir cada columna de la matriz es un vector de la base, se puede demostrar que T es otorgonal, y la transformación lineal definida por T se llama ortogonal tambien.


Si escribimos

                    |b11  b21  bn1|        |<v,b1>|

v*T = (v1,v2...vn)* |b12  b22  bn2|     =  |<v,b2>|

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

                    |b1n  b2n  bnn|        |<v,vn>|




o sea al mutiplicar un vector por la matriz T, que resulta de poner en cada columna un vector de la base ortonormal, se obtienen las coordenadas de v en la base B (Eso se suele notar como (v)B = v*T )

Nota: La coordenada i-esima es la proyección del vector sobre el vector i-esimo de la base.

Upgrade 2010:

En los graficos por computadora, es usual tener que transformar un punto del espacio real 3d (usualmente llamado world space) , al espacio del observador (usualmente llamado view space)
La función D3DXMatrixLookAtLH del DirectX calcula una matriz de vista, que luego se aplica a los puntos para transformarlos de world space al view space, en un proceso que se llama View Transform, en la jerga del DirectX

Esa matriz llamada View, esta compuesta por los vectores de una base ortonormal definidos por la direccion a la cual el observador esta mirando. El View Transform del pipeline del DirectX no es otra cosa que aplicar la expansion de Fourier en esa base a cada punto o vértice.
Usualmente se trabaja en en R4, para poder aplicar las traslaciones como multiplicación de matrices en la misma transformación. (Pero eso es tema del proximo post)


D3DXMATRIXA16 matView;
// posicion del punto de vista
D3DXVECTOR3 vEyePt(100,50,200);
// posicion a la que estoy mirando
D3DXVECTOR3 vLookatPt(0,0,0);
// direccion arbitraria hacia arriba
D3DXVECTOR3 vUpVec(0,0,1);
// el DX usa esta funcion para generar la matriz T
D3DXMatrixLookAtLH( &matView, &vEyePt, &vLookatPt, &vUpVec );


// segun lo que vimos, necesitamos una base ortonormal, que vamos a
// construi asi :
D3DXVECTOR3 vX,vY,vZ;
// la direccion Z es la direccion a la que estamos mirando
vZ = vLookatPt - vEyePt;
// la direccion X tiene que ser perpendicular a la direccion Z
// (el cross producto A = BxC devuelve un vector perpendicular a B y C,
// para eso necesitamos el vector UP, para poder generar un vector perpendicular
D3DXVec3Cross(&vX,&vUpVec,&vZ);
// pero la direccion UP, era aproximada, la real la generamos otra vez con
// el cross producto, Y es un vector perpendicular tanto a X, como Z
D3DXVec3Cross(&vY,&vZ,&vX);
// Hasta aca tenemos vX,vY,vZ todos perpendiculares entre si.
// nota: 3 vectores ortogonales son una base de R3
// Normalizamos
D3DXVec3Normalize(&vX,&vX);
D3DXVec3Normalize(&vY,&vY);
D3DXVec3Normalize(&vZ,&vZ);
// Ahora
// B = {vX,vY,vZ} = Base ortonormal de R3
// La matriz T = [vX vY vZ] es un aplicacion ortogonal

// Verificacion las columnas de T son los vectores de la Base
D3DXMATRIX T;
// vX vY vZ vW
T._11 = vX.x , T._12 = vY.x, T._13 = vZ.x, T._14=0;
T._21 = vX.y , T._22 = vY.y, T._23 = vZ.y, T._24=0;
T._31 = vX.z , T._32 = vY.z, T._33 = vZ.z, T._34=0;
T._41 = 0 , T._42 = 0, T._43 = 0, T._44=1;

// Hasta aca tenemos la parte de la matriz 3x3. Lo que falta para completar
// la matriz del DirectX que es de 4x4 es la translacion
T._41 = -D3DXVec3Dot(&vEyePt,&vX);
T._42 = -D3DXVec3Dot(&vEyePt,&vY);
T._43 = -D3DXVec3Dot(&vEyePt,&vZ);

// esto tambien se prodria haber hecho asi:
// | 1 0 0 0 |
// | 0 1 0 0 | x T
// | 0 0 1 0 |
// | -tx -ty -tz 1 |

// Trans = matriz de traslacion para origen al punto de vista.
// Entonces primero aplico la traslacion y luego el cambio de coordenadas
// pp dicho

//D3DXMATRIX Tras;
//D3DXMatrixIdentity(&Tras);
//Tras._41 = -vEyePt.x , Tras._42 = -vEyePt.y, Tras._43 = -vEyePt.z;
//D3DXMatrixMultiply(&T,&Tras,&T);

// La matriz T y la matView son iguales


Como decía anteriormente, la gracia de partir de una definición tan abastracta del producto interno radica en aplicarla a espacios no convencionales. Un buen ejemplo es el espacios de las funciones. Se puede demostrar facilmente que las funciones con las operaciones de suma y producto triviales (sumar funciones y multplicar una función por un numero), son un espacio vectorial.
Si pensamos una funcion como un vector de infinitas componentes, donde cada componente es el valor de la funcion en un punto, resulta intuitivo, que la sumatoria que define el producto interno habitual entre 2 vectores de dimensiones finitas, se transforma en una integral, lo que motiva a definir el producto interno entre 2 funciones como la integral del producto sobre el intervalo.
b
= integral f(x)g(x) dx
a

se puede demostrar que el conjunto de las funciones continuas acotadas en el intervalo [a,b], es un esp. vectorial y asi definido es un producto interno.

Con algunas condiciones más, las series de senos y cosenos típica de la forma de Fourier se basan en expresar una vector como una combinación lineal de elementos de la base, y en las propiedades de ortogonalidad

Nota: más tecnicamente se requiere un espacio de Hilbert, es un poco mas que un espacio vectorial.
El producto interno es una integral pero con algunas constantes para que cumpla las condiciones del producto interno, por ejemplo definida entre –pi y pi, la funcion tiene que ser cuadrado integrable etc etc) todas cosas que no vienen al caso en este contexto.

Como la base es infinita, si uno toma solo un cantidad finita de funciones (o elementos de la base), la serie se transforma en una suma, y la función resultante es una “proyección” de la función original en el sub-espacio generado por la base finita, y como habiamos dicho es la mejor aproximación que se puede obtener de esa función en dicha base. Es decir podemos ver a la suma de Fourier como una proyección en un espacio de dimension n, de la funcion real que vive en un espacio de dimension infinita.

Existen otras bases de espacios de funciones, como los polinomios ortogonales y que tienen propiedades similares a las de Fourier, usualmente se ven en la materia optativa de ecuaciones diferenciales de 3ero.

martes, 3 de agosto de 2010

Z-BUFFER Vectorial

El z-buffer es un invento magnífico, especialmente para todos aquellos que alguna vez tuvimos que resolver "a mano" el problema de las caras ocultas, para figuras arbitrarias. Cuando se trata de una superficie cerrada, se pueden ocultar facilmente las caras invisibles si se sabe de antemano en que orden están los vertices(horario o antihorario,pero todos iguales) . Si la proyección mantiene el orden la cara es visible (y si no esta oculta). (ver. Back-Face Culling). Pero para un conjunto arbitrario de caras, que es el caso más común en un escenario, el tema no es nada facil. Una de las formas de hacerlo es dibujando las objetos en orden z: primero los que estan mas alejados del pto de vista. Este algoritmo se llama del pintor. Para poder implementarlo es necesario una estructura especial con los caras, usualmente se utiliza en BSP (busquen binary space partition en el google, que hay cientos de ejemplos). En sintesis: es un despelote.

Todo esto se resuelve sencillamente con el directX activando el zbuffer, y con una resolucion de 24bits para la profundidad, es mas que suficiente para cualquier escenario tipico.

Ahora..si alguna vez intentaron trabajar con z-buffer y transparencias al mismo tiempo, la cosa se vuelve a complicar. Si bien en el directx permite configurar el alpha blending muy facilmente, el z-buffer asi como esta planteado no es compatible con el concepto de transparencia.
Una manera sencilla (y físicamente incorrecta) de resolverlo es dibujar primero todos los objetos opacos con el zbuffer activado, y luego desactivar el z-buffer y dibujar los objetos transparentes. Sin embargo, (y se puede demostrar facilmente) las ecuaciones que definen la transparencia NO SON CONMUTATIVAS, (salvo en el caso improbable que todos los objetos tengan el mismo color). En el caso general, la luz tiene que pasar por todos los objetos transparentes en el orden correcto (desde el más lejano hasta el más cercano) antes de alcanzar el punto de vista. Esto requiere que tengamos que dibujar los objetos transparentes en el orden z, por ejemplo usando una estructura de BSP, justamente lo que habíamos logramos evitar con el zbuffer!! . (El zbuffer surgió para evitar calcular estructuras como el BSP).

La solución que presento a continuacion es una prueba de concepto para extender el zbuffer. El zbuffer "normal" (y el color buffer) se puede pensar como "escalar" en el sentido que cada elemento del mismo contiene un valor único, que representa la profundidad, usualmente el valor de z que corresponde al punto más cercano al near plane.
Lo que necesitamos para implementar un alpha blending correcto, es que el zbuffer tenga todos los puntos al mismo tiempo, para poder ordenarlos por profundidad y asi generar la transparencia. Es decir un z-buffer "vectorial" y de ahi el titulo del post. La idea es simular este comportamiento "vectorial" de la siguiente manera:

1- Usar una textura auxiliares (con formato F32R) que va a representar el zbuffer, la vamos a llamar DepthMap.
2- Otra textura auxiliar (con formato RGB) va a representar el color buffer.
3- En el primer paso se dibuja toda la escena con el zbuffer activado, con un shader especial, para que grabe la profundidad en el depth map y el color en el colorbuffer. Este paso puede requerir 2 pasadas y 4 texturas, dependiendo de la implementación, pero esto no viene al caso. (Por ejemplo si no se usan multiples render targets). Luego de este paso, cada punto del depthMap y del ColorBufer contiene la profundidad y color respectivamente, del punto MAS cercano al near plane.
4- En el segundo paso, se vuelve a dibujar la escena, tambien con el z-buffer activado, pero con un pixel shader especial que hace uso del depthmap previamente generado y de la posibilidad que se introduce en el shader 3.0 de la semantica DEPTH. Esta permite que el pixel shader genere una profundidad como salida. (Siempre me pregunte para que podria llegar a servir eso...). El ps pregunta si la profundidad del pixel actual, es menor o igual a la que esta previamente grabada en el DepthMap. Si es asi, eso indica que ese punto ya fue "procesado", su color esta en el ColorBuffer y su profundidad en el DepthMap, Aprovechando que el ps puede cambiar la profundidad y que este punto no lo queremos volver a procesar, el ps asigna la profundidad en -1, haciendo al pto invisible. Y permitiendo que el directx ubique otro punto en dicho lugar. Si la profundidad por el contrario es mayor, se trata de un pto diferente, asi que devuelve la misma profundidad, para que el ztest sea el standard, y el color, combina el color actual con el que esta grabado en el colorbuffer, es decir hace un alpha blending a manopla.
Luego de este punto, en el dephmap esta grabado el SEGUNDO punto mas cercano al near plane, y en color buffer, la combinacion de ambos.
Este paso hay que volver a repetirlo una cierta cantidad de veces, hasta que finalmente no hay mas cambios en el depthmap y en el color buffer. En cada paso avanzamos "un lugar" hacia el farplane.

Para poder verlo graficamente, modifique el algoritmo para que borre (Clear) el colorbuffer en cada paso, asi se puede ver que puntos son visibles en cada paso, como el plano va avanzando desde el pto mas cercano hasta los mas lejamos.



Uniendo todos los colorBuffers:




Este es el PS que dibuja la escena en el ColorBuffer

void RenderScenePS( VS_OUTPUT In,float2 vpos:VPOS,out float4 Color:COLOR0,out float Depth:DEPTH)

{

   float prof = In.depth.x / In.depth.y; // prof del punto actual

   // pos = pos en el colorbuffer (0,1)

   float2 pos = float2(vpos.x/DX,vpos.y/DY);

   float prof_ant = tex2D(g_samDepthMap, pos).r;

   Color = 0;

   Depth = -1;

      if(prof>prof_ant+EPSILON )

      {

         // punto visible

         float4 scr_color = tex2D(g_samColorBuffer, pos);

         float4 dest_color = tex2D(MeshTextureSampler, In.TextureUV) * In.Diffuse;

         Color = scr_color + dest_color*kt;

      }

}




El ps que genera el depth map es muy similar, solo que en lugar del color devuelve la profundidad. Si se usaran MRT no haria falta 2 ps separados.



El algoritmo tiene 3 problemas importantes:
- El orden en que se obtienen los distintos mapas es inverso, primero el que esta mas cerca y por ultimo el que esta mas lejos. Eso es asi aproposito, primero probe de hacerlo al reves y me encontre con otro problema peor que no pude resolver. El problema que esten inversos se resuelve con un poco de algebra, al aplicar la ecuacion de alpha blending en cada paso, pero basicamente, el primer paso tiene mas peso que el segundo, y asi sucesivamente. Los ultimos pasos son los que menos contribuyen a la imagen final.
- No se sabe bien cuantos pasos hay que hacer, si bien "teoricamente" el algortimo se detiene cuando no hay cambios entre un paso y el siguiente, es muy costoso verificar eso en una textura de 1000 x 700. Por ese motivo decidi que el orden sea desde el near plane hacia el far plane, asi la cantidad de pasos la puedo dejar fija, y listo, total como los ultimos pasos son los que menos contribuyen, no seria tan grave perderlos. Si lo hiciera al reves, ahi si estaria en problemas...
- Hay serios problemas de redondeo, de momento hay que preguntar siempre con algun EPSILON puesto muy a manopla, con lo cual estoy trabajando para resolver estos problemas de inestabilidad del algoritmo.