Intersectando rayo con Cilindro

[ZüNdFoLGe]

Estoy tratando de implementar un raytracer, por ahora iba bastante bien hasta que me choque con el cilindro (hasta ahora llevo bien las demas primitivas: plano, esfera, cubo, triangulo).

El problema es que no doy pie con bola para que se vea el cilindro en pantalla (que se vea bien), con el codigo que muestro a continuacion el cilindro se ve punteado, sin color y con la sombra mal calculada (todos estos rollos son consecuencia de que hay un fallo en el calculo de la interseccion entre rayo y cilindro, que es con lo que estoy liado) . Tambien se esta renderizando como cilindro infinito (sin tapas) y en realidad lo estoy acotando (o lo estoy haciendo mal).

Me he guiado de este  tutorial (hasta ahora me ha sido muy util), cuyo autor creo que merodea estos foros (Juan Mellado):
http://www.inmensia.com/articulos/raytracing/cilindroycono.html?pag=1

Se que puede resultar bastante tedioso ponerse a ver un codigo de raytracing, pero si alguien por casualidad lo mira y puede echar una luz estare agradecido.

De la misma forma si alguien sabe de algun link con algun ejemplo sobre  interseccion entre rayo y cilindro se agradece que lo comparta (a pesar que ya he visto varios ejemplos y sigo sin dar pie con bola).

Código [Seleccionar] Expand


int Cylinder::IntersectRay( Ray& a_Ray,  float& a_Dist )
{
int retval = MISS;

float a = a_Ray.GetDirection().x * a_Ray.GetDirection().x + a_Ray.GetDirection().z * a_Ray.GetDirection().z;

float b = 2 * (a_Ray.GetOrigin().x * a_Ray.GetDirection().x + a_Ray.GetOrigin().z * a_Ray.GetDirection().z);

float c = (a_Ray.GetOrigin().x * a_Ray.GetOrigin().x) + ((a_Ray.GetOrigin().z * a_Ray.GetOrigin().z) ) - m_SqRadius;

        // calculo el determinante
float det = b*b - 4*a*c;

       
if (a != 0  &&  det >= 0) {  // si hay puntos de corte los calculo

float i1 = (-b - sqrtf(det))/(2*a);
float i2 = (-b + sqrtf(det))/(2*a);

if (det == 0) {
a_Dist = i1;
return HIT;
}

if (i1<0 && i2<0) return MISS;

                // obtengo el menor (positivo) de entre los 2 valores
float minor = getMinor2F(i1, i2);

// llamo "I" al punto de interseccion, calculo sus componentes a partir de la ecuacion del rayo
vector3 I = vector3(a_Ray.GetOrigin().x + minor* a_Ray.GetDirection().x,
a_Ray.GetOrigin().y + minor* a_Ray.GetDirection().y,
a_Ray.GetOrigin().z + minor* a_Ray.GetDirection().z
                   );

if (minor < a_Dist) {
// verifico las intersecciones en las "tapas" del cilindro
if ((I.x * I.x + I.z * I.z) <= 1 && (I.x * I.x + I.z * I.z) >= 0) {
a_Dist = minor;
retval = HIT;
}

}
}

return retval;
}


[nostromo]

Hola.

Busca en amazon http://www.amazon.com/Geometric-Computer-Graphics-Morgan-Kaufmann/dp/1558605940
Como sabrás, a veces puedes ver dentro del libro con "Search Inside". Busca "linear component cylinder" e intenta ir al apartado 11.3.4 pag. 508
Ahí te viene la explicación de un método de intersección linea-cilindro.

Te pongo un breve resumen del método que propone:

El método se basa en poner el cilindro en una forma estandar, a saber:
- Debes calcular una transformacion M para poner el cilindro con
 1. La base del cilindro asentada en el plano x,y
 2. El centro de la base alineado con el centro de coordenadas (0,0,0)
 3. El cilindro crece con altura h alineado a z
 4. El cilindro tiene radio r

Ahora, para calcular que un punto esta en el cilindro (sin contar las tapas):
Se tiene que cumplir (A)         x^2 + y^2 = r^2    y    0 <=z <= h

El rayo L(t)=P+t·D  tienes que transformarlo a las nuevas coordenadas del cilindro.... multiplicando por  M  

substituyes las componentes de L(t) en (A) y te sale un ecuación de segundo grado con la incognita t . Con la t calculada sacas z y compruebas que este en el rango 0 <=z <= h , (para que no este en el cilindro infinito).
Si el punto calculado esta en el cilindro se le aplica la transformación con la inversa de M para conseguir el punto en las coordenadas del mundo.

Para el calculo de la intersección con las tapas: calculas la intersección con el plano de la tapa y compruebas que  el punto (x,y)  cumple x^2 + y^2 <= r^2     Si cumple se hace lo mismo, se le aplica la transformación con la inversa de M y ya tienes el punto (x,y,z)

nota: En el libro se define M como la tranformación de las coordenadas de la  forma estandar a las coordenadas del mundo. Yo lo he puesto al revés, pero sigue siendo valido el resumen :)


Saludos

[Juan Mellado]

¡Que emoción, me han citado!  :o

Prueba a cambiar el m_SqRadius del final de esta línea:

Código [Seleccionar] Expand


float c = (a_Ray.GetOrigin().x * a_Ray.GetOrigin().x) + ((a_Ray.GetOrigin().z * a_Ray.GetOrigin().z) ) - m_SqRadius;


Por un 1:

Código [Seleccionar] Expand


float c = (a_Ray.GetOrigin().x * a_Ray.GetOrigin().x) + ((a_Ray.GetOrigin().z * a_Ray.GetOrigin().z) ) - 1;


Si dices que no tienes problema con planos, esferas y demás, voy a suponer que la transformación de unas coordenadas a otras lo tienes superado. En caso contrario puede ser más problemático.

Si quieres algún código de ejemplo puedes mirarte el mítico POV-Ray:
http://www.povray.org/redirect/www.povray.org/ftp/pub/povray/Official/Windows/povwin_s.zip
En la clase cones.cpp tienes las fórmulas para cilindro y cono. Pero ten en cuenta que la nomenclatura y los ejes de coordenadas son distintos a los mios (creo que Y es Z).

Lo más sencillo para estos casos es depurar con un ejemplo para un único rayo que sepas donde intersecta, en vez de lanzar toda la escena y ver que sale. Son números muy sencillos que puedes hacer a mano y comparar  depurando paso a paso con los valores de tu programa.

Una última cosa. Procura no comparar con cero (como en a != 0 por ejemplo). Por temas de precisión, lo mejor es utilizar un "epsilon", una cantidad muy pequeña que sea la base para la precisión con la que estás calculando (0.000001 por ejemplo, o algo similar)

Espero lo soluciones, los raytracers son muy majos.

[ZüNdFoLGe]

Primero que nada muchas gracias a ambos por las respuestas. Saque provecho de ambas y he avanzado, aunque todavia no se logra ver bien el cilindro.

De acuerdo al articulo en amazon que ha mencionado nostromo, el codigo lo tengo exactamente igual, excepto que el algoritmo de amazon multiplica por una matriz que representa  una transformacion inversa, la cual desconozco como esta compuesta. Esto lo hace antes de realizar los calculos:

Código [Seleccionar] Expand


Matrix4x4 transform, invTransform;
transform = cylinder.getTransformMatrix();
invTransform = transform.Inverse();
rayAuxiliar.direction = ray.direction * invTransform;
rayAuxiliar.origin = ray.origin * invTransform;



Seguramente mi error sea conceptual y precisamente este en esa parte de la transformacion de coordenadas, porque el resto del codigo esta igual que en el algoritmo antedicho.


Este es mi codigo:

Código [Seleccionar] Expand


int Cylinder::Intersect( Ray& a_Ray, float& a_Dist )
{

// origen del rayo, yo para "trasladarlo" al sist. de coord. del cilindro (que es donde posiciono al cilindro) le resto el centro del cilindro,
       //el centro es el punto (0, 0, 0), y obviamente la operacion no surte ningun efecto.
       // Aqui debe ser donde tengo el error pero no me doy cuenta cual es la transformacion que debo aplicarle al origen del rayo

       vector3 posRay = a_Ray.GetOrigin() - m_Centre;

// dir rayo, a la direccion del rayo no le estoy aplicando ninguna transformacion, la dejo como esta, pero sospecho que esto esta mal

vector3 dirRay = a_Ray.GetDirection();

       // de aqui para abajo esta todo igual al algoritmo de amazon que menciono nostromo

int retval = MISS;
float a = dirRay.x * dirRay.x + dirRay.z * dirRay.z;

float b =  (posRay.x * dirRay.x +
posRay.z * dirRay.z);

float c = (posRay.x * posRay.x) +
(posRay.z * posRay.z)  - 1;

float det = b*b - a*c;

if (a > EPSILON  &&  det >= EPSILON) {
float d = sqrt(det);
float i1 = (-b - d)/a;
float i2 = (-b + d)/a;

vector3 t0 = a_Ray.GetOrigin() + i1 * a_Ray.GetDirection();
// no valido
if (t0.y < EPSILON || t0.y > m_height) {
valid[0] = false;

} else {
valid[0] = true;
}

vector3 t1 = a_Ray.GetOrigin() + i2 * a_Ray.GetDirection();
if (t1.y < EPSILON || t1.y > m_height) {
valid[1] = false;
} else {
valid[1] = true;
}

num[0] = i1;
num[1] = i2;
float minor;
if (valid[0] && valid[1]) minor = getMinValue(num[0], num[1]);
else {
for(int i=0; i<2; i++) {
if (valid[i]) minor = num[i];
}
}

if (minor < a_Dist) {
retval = HIT;
a_Dist = minor;
}
}

return retval;
}


El cilindro lo estoy dibujando con estos datos:
Radio = 1
Centro = (0, 0, 0)
Altura (h) = 6

Aun sigo luchando y estoy leyendo este articulo de Juan Mellado http://www.inmensia.com/articulos/raytracing/interseccion.html?pag=2 , pero he estado muchas horas con esto y seguramente mi error sea conceptual.

[nostromo]

Hola,

Te falta multiplicar por 2:
float b =  2 * (posRay.x * dirRay.x + posRay.z * dirRay.z);

Se te olvido multiplicar por 4:
float det = b*b - 4*a*c;

Se te olvido dividir por 2a:
                float i1 = (-b - d)/2*a;
      float i2 = (-b + d)/2*a;

Creo que te vendria bien hacer, para que tu lo visualizes mejor, un esquema con las coordenadas y los ejes de tu mundo.
Viendo esta linea de codigo:   "if (t0.y < EPSILON || t0.y > m_height) { ......."  entiendo que tu cilindro "crece en altura" en el eje y.
El caso es que las formulas estan hechas asumiendo que el cilindro crece en el eje z y tenga una base(tapa) en el plano x,y.

Yo de ti primero probaria sin transformar nada. Hasta que te salga algo coherente. Y cuando te salga... entonces metete con la transformacion

[ZüNdFoLGe]

Finalmente he podido renderizar el cilindro, me he bajado el codigo del POV-ray para fijarme como intersecta el rayo con las cuadricas y de ahi adaptarlo a lo que ya tenia. Aqui va la imagen, por defecto el observador esta en el punto (0, 0, 0):



Ahora estoy con un nuevo problemita, y es que cuando muevo la camara la direccion del rayo no esta saliendo de "el ojo del observador", a pesar de que a mi entender estoy actualizando bien la direccion a medida que muevo la camara. Asi es como se ve la imagen cuando coloco al observador en el punto (-400, 0, 0):



Asi es como estoy seteando a la direccion del rayo y la del observador:

Código [Seleccionar] Expand


    for(int y = 0; y < IMAGE_Y; y++) {
        for(int x = 0; x < IMAGE_X; x++) {
               
            position observador,  direccion;

            observador.x = -400;
            observador.y = 0;
            observador.z = 0;

            direccion.x = (WINDOWS_X/IMAGE_X)*x - WINDOWS_X/2;
    direccion.y = (WINDOWS_Y/IMAGE_Y)*y - WINDOWS_Y/2;
           
            direccion.z = WINDOWS_Z;

     color col = traza_RayTracing(observador, direccion, 1,-1,objects,lights,count,countLights);


Aparentemente viendo como queda la imagen que les mostre (que la esfera queda como un huevo, como si fuera un elipsoide) cuando yo cambio la posicion del observador (lo que hago es cambiar observador.x, observador.y) la direccion del rayo sigue partiendo del 0, 0, 0... creo que ese es el unico motivo para que las imagenes se vean "ahuevadas" como muestra la figura.

Si pruebo restar direccion - observador para obtener la direccion "actualizada" no logro ver los objetos, la imagen se va de foco. Evidentemente algo estoy haciendo mal. Si alguien detecta que estoy cometiendo algun error en el calculo de la direccion en el codigo que adjunte agradezco una mano.

Gracias de antemano,