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.