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 :