Derivar la Matriz de Proyeccion

Hay que eliminar una dimension.

Para poder dibujar en la pantalla de una computadora un mundo en 3d es necesario de alguna manera transformar los puntos del espacio real en 3 dimensiones, al llamado espacio de la pantalla, que logicamente tiene 2 dimensiones.
Este proceso de eliminar una dimension en algebra se llama proyeccion.
Por ejemplo, si defino

P:R3-->R2 / P(x,y,z) = (x,y)

es una proyeccion.

El problema con esta proyeccion es que no brinda demasiada informacion acerca del mundo 3d real. Eso es debido a que la dimension z simplemente desaparecio. Seria bueno de alguna manera que parte de la informacion que esta en esa coordenada aparezca en la funcion.

Proyeccion en perspectiva.

La perspectiva es un efecto visual, por el cual los objetos que estan cerca se ven mas grandes que los que estan lejos.

Para entender como funciona la perspectiva, es util el siguiente esquema:


En el esquema el near plane corresponde a la pantalla de la computadora, que esta a una cierta distancia "d" del ojo. El origen de coordenadas (0,0,0) esta en el ojo, que tambien podria ser un camara, por lo cual se lo llama camera space usualmente.
En la imagen solo esta la coordenada y, (y la z que es la distancia al ojo), para simplificar el esquema, pero la coordenada x funciona de la misma manera.

Volviendo a la perspectiva, si quisieramos que lo que esta mas lejos se viera mas pequeño, podriamos dividir por z y listo:

P:R3-->R2 / P(x,y,z) = (x/z,y/z)

Esta nueva proyeccion si bien esta bastante cerca de lo que queremos lograr, todavia tiene algunos problemas: en la medida que la z se hace cada vez mas pequeña, el valor de la proyeccion tiende a infinito. Por otra parte no podemos permitir que la z sea infinitamente grande, debido a la precision limitada de la computadora.

Clipping.

Entonces, lo primero que tenemos que hacer es acotar los valores de z entre ciertos limites que aseguren que exista la proyeccion y que se pueda mapear con una determinada precision.
Clipping es el proceso por el cual se elimina o se descartan los valores que estan fuera de cierto rango. En nuestro caso, la variable z tiene que ir entre el near plane y el far plane. (Si bien el far plane lo podriamos ubicar arbitrariamente lejos del near plane, dentro de los parametros logicos)

Volviendo a la imagen, el punto p = (x,y,z) que esta en el mundo real, y expresado en camera space (o coordenadas del ojo) , y lo queremos proyectar en la plantalla de la computadora, que es el near plane, y en el esquema esta representado por ps = (xs,ys) (que seria la coordenada en screen space).

Por triangulos semejantes, se tiene que
y/z = ys/d

con lo cual
ys = d * y/z

Ahora, queremos "normalizar" estos valores, para que vayan entre -1 y 1, de esta forma se pueden adaptar a cualquier resolusion de pantalla (800x600, 1024x768 etc)

Entonces:

H = Alto de la pantalla, y d = znear (como se puede ver en el esquema)
ys = y * (2 * znear ) / ( z * H)

Notar que la formula sigue siendo ys = y/z * K

de la misma forma
xs = x * (2 * znear ) / ( z * W)
donde W = ancho de la pantalla

Como cualquier cosa que se vaya de pantalla no se va a ver,
resulta que -1<=xs<=1 ; -1<=ys<=1

Que hacemos con la z?

La z se comporta diferente de las otras 2. Primero porque como no se puede ver nada que este detras de la pantalla, y por convencion zs va entre cero y 1, y no entre -1 y 1 como xs e ys.
El otro motivo tiene que ver con las:

Coordenadas Homogeneas.

Como todas estas ecuaciones tienen que ir a parar a una matriz de 4x4 para poder utilizar el hardware de la placa de video, que tiene implementadas multiplicaciones de matrices. Y como tenemos que dividir tanto la x como la y, por el valor z (constante mas, constante menos), aca entran en juego las coordenadas homogeneas.

En las coordenadas homogeneas todo punto en espacio real (x,y,z) esta definido por 4 coordenadas

(x * w ,y * w ,z * w , w) con w !=0

dichas coordenadas se llaman homogeneas y transforman el espacio R3 en el espacio proyectivo, pero eso es otro tema (de topologia)

Para obtener las coordenadas reales hay que dividir por w, es decir si (X,Y,Z,W) son homogeneas, para obtener x real, se hace x = X/W y = Y/W z = Z/W
En concreto es un artificio algebraico para poder introducir una division en una multiplicacion de matrices.

Nota: agregar una dimension adicional permite introducir traslaciones junto con rotaciones y escalados en la misma matriz.

Poniendo X e Y en la matriz de Proyeccion.

Antes de seguir con Z, vamos a ver que forma va tomando la matriz de proyeccion:

            | 2*Zn/W         0         0     0 | 
(x,y,z,1) x | 0             2*Zn/H     0     0 | 
            | 0             0         A?     1 | 
            | 0             0         B?     0 | 

Vemos que hasta aca se verifican:

X = 2*x*Zn/W
Y = 2*y*Zn/H
Z = Az + B
W = z

Haciendo la division por W, nos queda
xs = x * (2 * Zn ) / ( z * W)
ys = y * (2 * Zn ) / ( z * H)

Tal como queriamos, hasta aca estamos bien.
Ahora, para la Z podemos tocar las constantes A y/o B en la matriz, sin romper nada de lo que estaba.

Z = Az + B

Y luego de dividir por W=z, nos quedaria:

Zs = A + B/z

Pero sabemos que si
 z = Zn implica que  Zs = 0 y si
 z = Zf  impplica que Zs = 1, luego

0 = A +B/Zn
1 = A +B/Zf

Se trata de 2 ecuaciones con 2 incognitas (A y B), que despejando se resuelven asi:

A = Zf / (Zf-Zn)
B = Zf*Zn / (Zn-Zf)

Entonces, la Zs nos quedaria asi:

Zs = Zf/(Zf-Zn) + Zf*Zn / (z * (Zn-Zf))

Poniendo todo junto.

        | 2*Zn/W         0         0             0 | 
MProj = | 0             2*Zn/H     0             0 | 
        | 0             0         Zf/(Zf-Zn)     1 | 
        | 0             0         Zf*Zn/(Zn-Zf)  0 | 

Semi-Interactive Rendering

Entre el real time rendering dibujando a mas de 30 fps y las imagenes casi fotograficas generadas con los raytracer modernos, que llegan a tardar mas de 30min en una maquina de ultima generacion, se habre una brecha no muy explotada de aplicaciones que precisan crear graficas de la mayor calidad posible en un tiempo razonable que van desde uno pocos segundos hasta algunos minutos. Sin las limitaciones del realtime, se puede aprovechar la gpu para lograr efectos de sombras, reflexiones, perceptualmente correctas.

GPU Compute: Cuda del Subdesarrollo.

La semana pasada fui a una charla de Programación de Procesadores Gráficos Nvidia con Cuda que dieron en el dpto de electronica de la UTN. Estuvo espectacular porque además de que el flaco que dió la charla la tenía reclara, mostró ejemplos concretos de programas y comparaciones de algoritmos programados en cpu y gpu. Y por último como optimizar código para una arquitectura en particular. Es fantástico ir a una presentación donde aparece código fuente en las diapositivas!

A simple vista me pareció bastante similar a los "hacks" que estamos habituados a hacer para meter cálculo en el Pixel shader. Por lo cual, me decidí a subir al blog esta prueba de concepto de cómputo usando la GPU, sin cuda, en directX 9: Cuda del Subdesarrollo.

El pequeño ejemplo simula un conjunto de bolas que ruedan sin deslizar en una superficie teniendo en cuenta sólo el efecto de la gravedad, un rozamiento y sin tener en cuenta las colisiones. Bajo esas hipótesis, el movimiento de cada bola es independiente del resto, por lo cual, es sumamente paralelizable.

La idea central es que el código que calcula y actualiza la posición de una bola específica se ejecuta en paralelo en cada pixel shader. Para ello tenemos varios problemas:

1- Como hacer que se ejecuta el PS una vez para cada bola.
2- Como sabe el PS que bola especifica tiene que procesar
3- Como usar la memoria para almacenar datos intermedios y resultados.


1- Como hacer que se ejecuta el PS una vez para cada bola.

El primer punto es muy sencillo, siendo que la placa lo único que entiende son triángulos, hay que dibujar 2 de ellos que formen un "quad", del tamaño adecuado, en general de un tamaño N x N = N^2 bolas. Por ejemplo, un quads de 8 x 8, podria generar un cuadrado compuesto de 64 pixeles. Y en ese caso llamaria al PS 64 veces.

CustomVertex.PositionTextured[] vertices = new CustomVertex.PositionTextured[]
{
new CustomVertex.PositionTextured( -1, 1, 1, 0,0),
new CustomVertex.PositionTextured(1, 1, 1, 1,0),
new CustomVertex.PositionTextured(-1, -1, 1, 0,1),
new CustomVertex.PositionTextured(1,-1, 1, 1,1)
};

La transformación tiene que ser tal, que se aplique un pixel por unidad, es decir no tiene que transformar ningun punto.

device.Transform.Projection = Matrix.Identity;
device.Transform.World = Matrix.Identity;
device.Transform.View = Matrix.Identity;

2- Como sabe el PS que bola especifica tiene que procesar

Ya que el código del PS es el mismo para todas las bolas, dentro del mismo hay que poder identificar cual de las bolas estamos calculando. Una manera muy fácil es usar las coordenadas de textura. Como se puede ver el Vertex buffer está creado de tal manera que las coordenadas uv mapean uniformemente todos los puntos: (0,0) (1,0), (0,1) y (1,1).


3- Como usar la memoria para almacenar datos intermedios y resultados.

Aprovechando que cada pixel tiene un ID u,v de las coordenadas de texturas, toda la información adicional y demás resultados, se almacenará en texturas de NxN.
En el ejemplo se usan texturas para:

Velocidad:
g_pVelocidad = new Texture(d3dDevice, MAX_DS, MAX_DS, 1, Usage.RenderTarget,
Format.A32B32G32R32F, .Default);

Posición:
g_pPos = new Texture(d3dDevice, MAX_DS, MAX_DS, 1, Usage.RenderTarget,
Format.A32B32G32R32F, Pool.Default);

Normal:
g_pNormal = new Texture(d3dDevice, MAX_DS, MAX_DS, 1, Usage.RenderTarget,
Format.A32B32G32R32F, Pool.Default);

Tangente:
g_pBiNormal = new Texture(d3dDevice, MAX_DS, MAX_DS, 1, Usage.RenderTarget,
Format.A32B32G32R32F, Pool.Default);

Además de otras varias temporarias, y como no se puede leer y escribir sobre la misma textura al mismo tiempo, cada textura tiene su contrapartida Out:
g_pVelocidadOut, g_pPosOut, etc.

Y al final de cada paso, se intercambian las texturas, la "out" pasa a ser "in", y viceversa:

Texture aux = g_pVelocidad;
g_pVelocidad = g_pVelocidadOut;
g_pVelocidadOut = aux;


Dentro del contexto de un PS, la coordenada uv de textura permite indexar todos los valores necesarios, como la velocidad y la posición. Y como la salida está redirigida también a una textura de salida, dicho valor siempre se almacena en la misma posición.


Los resultados se almacenan aprovechando la salida del PS, que en lugar de ser un "color" representa en cada paso 4 valores float. Que pueden representar posicion, velocidad, normales etc. Como hace falta mas que 4 valores hay varias pasadas, y los valores float se van cableando en distintos canales de la textura.
Hay que notar que la coordenada UV coincide con la posicion en pantalla, esta todo alineado de tal manera que cada PS tiene asigando el mismo lugar en las texturas y en la salida. Por ejemplo el pixel de arriba a la izquierda tiene las coordenadas u=0, v=0. Y por lo tanto la posicion actual, la puede tomar de la textura de posiciones en u=0, v=0. Y la salida, tambien va a ir a parar la textura de posicionOut a la posicion u=0, v=0.

Ejemplo paso x paso:


En el Ps, la técnica tiene 3 pasadas, una para inicializar posición y velocidad, otra para calcular pp dicha, que es la que se va a ejecutar en cada render, y otra de preview que es para visualizar los resultados. Las 2 primeras son las importantes.

// Genera el mapa de velocidad
technique ComputeVel
{
    pass P0
    {          
        PixelShader  = compile ps_2_0 PSInitVel(); 
    }
    pass P1
    {          
        PixelShader  = compile ps_2_0 PSComputeVel(); 
    }
    
    pass P2
    {          
        PixelShader  = compile ps_2_0 PSPreview(); 
    }
    
}

El primer paso es inicializar las texturas de posición y velocidad. O sea lo que seria la posición y velocidad inicial de cada bola a partícula. Para ello voy a usar 2 render target al mismo tiempo (si la placa no se lo banca hay que hacer mas pasos y listo).

            // rt1 = velocidad
            pOldRT = device.GetRenderTarget(0);
            pSurf = g_pVelocidad.GetSurfaceLevel(0);
            device.SetRenderTarget(0, pSurf);
            // rt2 = posicion
            Surface pSurf2 = g_pPos.GetSurfaceLevel(0);
            device.SetRenderTarget(1, pSurf2);

            device.Clear(ClearFlags.Target | ClearFlags.ZBuffer, Color.Black, 1.0f, 0);
            Effect g_pEffect = mesh.Effect;
            device.BeginScene();
            g_pEffect.Technique = "ComputeVel";
            device.VertexFormat = CustomVertex.PositionTextured.Format;
            device.SetStreamSource(0, g_pVBV3D, 0);
            g_pEffect.Begin(FX.None);
            g_pEffect.BeginPass(0);
            device.DrawPrimitives(PrimitiveType.TriangleStrip, 0, 2);
            g_pEffect.EndPass();
            g_pEffect.End();
            device.EndScene();
            device.Present();

Definimos una estructura auxiliar, que permite devolver al mismo tiempo los valores de pos. y vel:

struct PS_OUTPUT
{
    float4 Velocity : COLOR0;
    float4 Position : COLOR1;
};

El PS pp dicho:

PS_OUTPUT PSInitVel( float2 Texcoord: TEXCOORD0)
{
    PS_OUTPUT Output;
    Output.Velocity = float4(0,0,0,0);
    
    float x0 =     0.2 + 0.6*(Texcoord.x-0.5)*currentScaleXZ*map_size;
    float z0 =  0.2 + 0.6*(Texcoord.y-0.5)*currentScaleXZ*map_size;
    Output.Position = float4(x0,z0,2000,0);    
    return Output;
}

Esta es una funcion auxiliar, que permite calcular la altura del heighmap en un determinado punto, usando un texture look up:

float CalcularAltura(float x,float y)
{
    float u = y/(currentScaleXZ*map_size) + 0.5;
    float v = x/(currentScaleXZ*map_size) + 0.5;
    float3 T0 = tex2D(heightMap,float2(u+map_desf,v+map_desf));
    float H0 = T0.r*0.299 + T0.g*0.587 + T0.b*0.114;
    return 255*H0*currentScaleY;
}

En cada paso de la simulación, que se llama cada vez que hacemos un OnRender(), se ejecuta el PS PSComputeVel, que calcula la nueva velocidad y actualiza la posicion y velocidad para el siguiente paso:

            // rt1= velocidad
            pOldRT = device.GetRenderTarget(0);
            pSurf = g_pVelocidadOut.GetSurfaceLevel(0);
            device.SetRenderTarget(0, pSurf);
            // rt2 = posicion
            Surface pSurf2 = g_pPosOut.GetSurfaceLevel(0);
            device.SetRenderTarget(1, pSurf2);

            device.Clear(ClearFlags.Target | ClearFlags.ZBuffer, Color.Black, 1.0f, 0);
            device.BeginScene();

            Effect g_pEffect = mesh.Effect;
            mesh.Effect.SetValue("elapsedTime", elapsedTime);
            g_pEffect.Technique = "ComputeVel";
            mesh.Effect.SetValue("g_pVelocidad", g_pVelocidad);
            mesh.Effect.SetValue("g_pPos", g_pPos);
            device.VertexFormat = CustomVertex.PositionTextured.Format;
            device.SetStreamSource(0, g_pVBV3D, 0);
            g_pEffect.Begin(FX.None);
            g_pEffect.BeginPass(1);
            device.DrawPrimitives(PrimitiveType.TriangleStrip, 0, 2);
            g_pEffect.EndPass();
            g_pEffect.End();
            device.EndScene();

            // swap de texturas
            Texture aux = g_pVelocidad;
            g_pVelocidad = g_pVelocidadOut;
            g_pVelocidadOut = aux;
            aux = g_pPos;
            g_pPos = g_pPosOut;
            g_pPosOut = aux;

PS_OUTPUT PSComputeVel( float2 Texcoord: TEXCOORD0)
{
    PS_OUTPUT Output;
    float4 Vel = tex2D( Velocidad, Texcoord );
    float vel_x = Vel.x;
    float vel_z = Vel.y;
    
    float4 Pos = tex2D( Posicion, Texcoord );
    float x0 = Pos.x;
    float z0 = Pos.y;
    float H = CalcularAltura(x0, z0);
    float ddx,ddz,alfa;
    if(Pos.z>H+10)
    {
        H = Pos.z - elapsedTime*500;
        ddx = ddz = 0;
        alfa = 0;
    }
    else
    {
        float dt = 20;
        ddx = (CalcularAltura(x0 + dt, z0) - H) / dt;
        ddz = (CalcularAltura(x0, z0 + dt) - H) / dt;
        vel_x -= elapsedTime * ddx * 150;
        vel_z -= elapsedTime * ddz * 150;
        
        // rozamiento
        vel_x *= 1-elapsedTime/25;
        vel_z *= 1-elapsedTime/25;

        x0 += elapsedTime * vel_x;
        z0 += elapsedTime * vel_z;
    
        // calculo el angulo rodado:
        float dx = x0 - Pos.x;
        float dz = z0 - Pos.y;
        alfa = sqrt(dx*dx + dz*dz) / Kp;
    }
    
    
    Output.Velocity = float4(vel_x,vel_z,ddx,ddz);
    Output.Position = float4(x0,z0,H,alfa);    
    return Output;
}

Aca lo importa es notar que Texcoord identifica una bola o partícula dentro de todas las que hay. Hay que pensar el que conjunto de particulas como una grilla o matriz de N x N, a la que se accede con la coordenadas de Texcoord, como el Matrix[i,j] de un programa en C.
Por ejemplo, para acceder a la posicion actual, se hace

float4 Pos = tex2D( Posicion, Texcoord );
float x0 = Pos.x;
float z0 = Pos.y;

La salida, va a ir a parar a la misma posición en la textura de salida. Esta es quiza una de las mayores desventajas con respecto al Cuda. Este hack no puede hacer que la salida del PS vaya a parar a otro lugar. Si quisiera modificar desde este PS, la posición de otra partícula, asi como esta hecho no se podría, porque la salida siempre va a ir a parar a la misma coordenada de textura que la entrada. (O a otra arbitraria, pero siempre a la misma). De momento no encontré forma de "cambiar" la posición del punto dentro del PS. Para hacer algo así deberiamos usar un VertexShader, que ahi si, "tocando" las coordenadas homogeneas, se puede hacer caer la salida en cualquier lado. Un aproach de ese estilo, lo use en el real time radiosity.

Dejo unas imágenes con 512 x 512 partículas y para la próxima van comparaciones con el mismo ejemplo en CPU:

Reflexiones planas usando el Stencil. (Espejos)

Como comentamos en otros post, las técnicas de environment map, (ya sea cube map o sphere maps) no se llevan bien con las reflexiones planas: los objetos del tipo espejos, o superficies planas brillantes, como pisos encerados, o mesadas, son muy sencibles a los artifacts de los enviroments maps. Para peor, el ojo humano esta muy acostumbrado a ver reflexiones en los espejos, y nota rapidamente cuando las mismas estan mal.

Uno de los métodos más usados para dibujar espejos consiste en re-dibujar la escena reflejada. Es decir dibujar todo otra vez, pero cambiando el look from y el look at, de tal manera de obtener la escena "reflejada".



El problema es que no queremos que el dibujo se "salga" de los limites del espejo o de la superficie en cuestion.
Muchas veces, esto se logra dibujando las cosas en un orden preciso: primero la escena reflejada, luego la escena real. En el post anterior vimos como el juego de tron usaba esa tecnica para la reflexion en el piso.

En el caso general esto no se puede hacer o no se puede asegurar, y la técnica mas comun para que el dibujo no se pase de los limites del espejo, es usar el stencil buffer, que funciona como un "máscara" que indica que pixeles se tienen que dibujar y cuales no.

Este metodo es el mas general, ya que por usar el stencil como mascara, soporta cualquier tipo de forma que tenga el "espejo": por ejemplo el caso de un piso de baldosas, donde solo algunas baldosas reflejan y otras no. O una mesada, donde la textura indica si hay o no reflexion. Esto gracias al que el stencil puede controlar que puntos se dibujan o no a nivel de pixel.

La técnica en si es muy sencilla, en verdad el desafio es reducir la cantidad de geometria que se renderiza en las reflexiones, con el uso de otras tecnicas de optimizacion que no vienen al caso.

Paso x paso seria algo asi: (nota: A veces conviene dibujar los espejos primero, aunque no es indispensable).

1- Antes de dibujar el espejo, hay que crear y preparar el stencil.

    // Guardo el zbuffer actuales
    LPDIRECT3DSURFACE9 pDSOld = NULL;
    GetDepthStencilSurface( &pDSOld);
    // Creo el stencil
    CreateDepthStencilSurface(d3dpp.BackBufferWidth,d3dpp.BackBufferHeight,D3DFMT_D24S8,D3DMULTISAMPLE_2_SAMPLES,0,TRUE,&g_pStencil,NULL);
    SetDepthStencilSurface( g_pStencil);
    Clear( 0, NULL, D3DCLEAR_ZBUFFERD3DCLEAR_STENCIL,0,1,0);

2- Hay que inicializar el stencil, poniendo "1" en los pixeles donde va el espejo, para ello hay que:

    // activar el stencil
    SetRenderState( D3DRS_STENCILENABLE,TRUE );

    // queremos que siempre escriba un 1
    SetRenderState( D3DRS_STENCILFUNC,D3DCMP_ALWAYS);
    SetRenderState( D3DRS_STENCILREF,1);
    SetRenderState( D3DRS_STENCILPASS,D3DSTENCILOP_REPLACE);

    // no queremos que escriba ningun color, no interesa el colorbuffer, solo queremos inicializar el stencil
    SetRenderState( D3DRS_COLORWRITEENABLE,0x00000000);
    SetRenderState( D3DRS_ZWRITEENABLE,FALSE);
    SetRenderState( D3DRS_ZENABLE, FALSE);
    // tampoco queremos texturas,
    SetTexture( 0, NULL);

    // dibujo el espejo pp dicho: (se supone que estas primitivas dibujan 2 triangulos que representan el espejo)
    DrawPrimitive( D3DPT_TRIANGLEFAN,0, 2);

    // luego de esto queda el zbuffer limpio, el color buffer limpio
    // y el stencil con 1 en los pixels que tienen espejo
    // dejo todo como estaba antes:
    SetRenderState( D3DRS_ZWRITEENABLE,TRUE);
    SetRenderState( D3DRS_ZENABLE, TRUE);
    SetRenderState( D3DRS_STENCILPASS,D3DSTENCILOP_KEEP);
    SetRenderState( D3DRS_COLORWRITEENABLE,0x0000000f);
    SetRenderState( D3DRS_STENCILFUNC,D3DCMP_LESSEQUAL);

3- Ahora tengo que dibujar la escena, como se veria reflejada, afortunadamente, hay una funcion de utilidad del directx que devuelve la matriz de reflexión.
Solo necesita el plano sobre el que está el espejo y hay una estructura en directX para definir un plano:

    // Preparo el lano del espejo Pn = normal al plano del espejo, plano_D = distancia al plano
    D3DXPLANE espejo;
    espejo.a = Pn.x;
    espejo.b = Pn.y;
    espejo.c = Pn.z;
    espejo.d = plano_D;

    // Preparo la matriz de reflexion
    // ------------------------------
    D3DXMatrixLookAtLH( &matView, &vEyePt, &vLookatPt, &vUpVec );
    D3DXMATRIX mReflectView;
    D3DXMatrixReflect(&mReflectView,&espejo);
    D3DXMatrixMultiply(&matView, &mReflectView, &matView);
    SetTransform( D3DTS_VIEW, &matView );

4- Ahora, esto es un poco mas sutil, solo hay que dibujar las cosas que estan delante del espejo!!.




Afortunadamente, tambien hay una forma de decirle eso al directx con lo que se llaman Clip Planes, solo hay que pasarle el plano, que en principio es el plano del espejo.

SetClipPlane(0,espejo);

El problema es que a veces queda al reves!! , y dibuja lo que esta atrás del espejo. En ese caso hay que multiplicar por -1. Para calcular si hay que multiplicar por -1 o no, hay que ver de que lado del plano queda el punto de vista, (el vEyePt), para ello primero hay que transformarlo al espacio de la matriz de vista reflejada, y usar esta formula que si es negativa o positiva o cero, nos indica de que lado del plano esta, o si esta justo en el plano.

ndist = Pn.x*P.x + Pn.y*P.y + Pn.z*P.z + plano_D;

Esta fórmula se interpreta asi:
Si ndist<0 el punto P esta de un lado del plano, Si ndist>0 el punto P esta del otro lado,
y si ndist es cero, esta justo sobre el plano.


Aplicando esta fórmula entonces:

    D3DXVECTOR4 pOut;
    D3DXVec3Transform(&pOut,&vEyePt,&mReflectView);
    LF.y = pOut.x;
    LF.x = pOut.y;
    LF.z = pOut.z;
    double K = ndist(TVector3d(LF.y,LF.x,LF.z))<0?1:-1;
    float plano[] = {Pn.x*K,Pn.y*K,Pn.z*K,plano_D*K};
    SetClipPlane(0,plano);
    SetRenderState( D3DRS_CLIPPLANEENABLE,TRUE);

5- Dibujamos toda la escena pp dicha. (Menos el espejo!! logicamente) y una vez finalizado, dejamos todo como estaba antes de empezar la reflexion:

    // restauro el stencil global de la escena
    SetDepthStencilSurface( pDSOld);
    // deshabilito el stencil
    SetRenderState( D3DRS_STENCILENABLE,FALSE);
    // restauro el render state 
    SetRenderState( D3DRS_ZENABLE, TRUE);
    SetRenderState( D3DRS_ZFUNC, D3DCMP_ALWAYS);

6- Ya "casi" estamos, ahora, un "snapshot" del color buffer, mostraría la escena dibujada sobre el espejo, y solamente sobre el espejo, ya que eso lo logramos con el stencil. Pero el zbuffer no esta correcto, ya que nosotros dibujamos en un depthbuffer auxiliar. Hay que tener en cuenta que el stencil buffer se almacena junto con el zbuffer, en lo que se llama DepthStencilBuffer.
Lo que necesitamos es que la profundidad grabada NO SEA la del dibujo, si no la del espejo. (Como se fuera opaco). El problema es que el espejo refleja los colores de un objeto pero el z tiene que seguir siendo el del espejo, no el del objeto reflejado. De esta forma cuando el resto de la escena se dibuje funcione correctamente el ztest.
Eso lo resolvemos simplemente dibujando el espejo, totalmente transparente, y con el ztest activado.


Por otra parte, el espejo se podria dibujar parcialmente transparente, inclusive con una textura para dar un efecto de superficie espejada, y no espejo al 100%. Es decir combinar ambas cosas, como una superficie pulida, como el piso de la siguiente escena:





Tambien, se pueden combinar varios espejos aplicando la misma tecnica varias veces




Como el stencil tambien se puede configurar para que sume, o reste cada vez que se dibuja algo, es posible, modificar esta misma técnica para que genere reflexiones "recursivas". Es decir que si mientras esta dibujando un espejo, se ve otro espejo, tiene que dibujar otra escena reflejada de la reflejada, y asi hasta una cierta profundidad en una especie de "arbol" de reflexiones. Para ello se saca provecho del "contador" que tiene el stencil, y cada vez que se dibuja un espejo, en lugar de escribir un "1", se le indica que "sume 1" al valor que tenia antes. De esta manera se van recortando los distintos espejos, de forma recursiva. Cuando termina de dibujar un espejo, para restaurar el stencil a como estaba antes, vuelve a dibujarlo pero le dice al stencil que "reste 1".
No me consta si esto es algo standard o no, porque si bien lo pude hacer funcionar (casi bien), no resultó práctico por la cantidad de veces que dibuja la escena, que crece en forma exponencial a la cantidad de espejos y de reflexiones.
Probablemente los motores que muestran muchas reflexiones recursivas esten usando una tecnica diferente.

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:

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 :

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.