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 :