Lyapunov 3D

English    Spanish

Today I'm going to comment a small program that it's not fully related to raytrace, but is related to 3D.

It has been some time since I have tried to represent the Markus-Lyapunov fractal. The 3D version extends the exponent computation one dimension. I have used CUDA and OpenGL to represent that.

Thanks to CUDA it's possible do computation using the GPU, so, if you have a good video card, the computation will be quite fast.

The variation of the algorithm is simple. The sequences contain one more letter, in this case C, that represents a point in the Z coordinate (or 'c'). Everything else is similar to the previous algorithm.

In order to do the renders I use two VBOs: One for storing positions and another one to store colors. The positions will be arranged in a volume of 256x256x256 elements. To compute the colors of each voxel I use two kernels. One computes the positions of each point and stores them in one VBO. The other use that positions to compute the color, that is the lyapunov exponent in that point. The points with an exponent greater than zero (or greater than a cut-off value) are moved to a region in order to not disturb. The points with a negative exponent (or lower than the cutoff value) will be drawn with a color proportional to the absolute value of the exponent.

Space region (0,4)-(0,4)-(0,4)

Space region (1,4)-(1,4)-(1,4)

Adjusting the cut-off value in order to remove some information about stability.

Readjusting the cut-off value to be nearer to zero.

The code for computing the Lyapunov 3D exponent for the sequence ABC is this:

__device__ float evalLyap(float3 p)
{
float sum=0;
float Xn=0.5f;


for (int i = 0; i < 3; i++) {
 float rn = i==2 ? p.x : i==1? p.y : p.z;
 Xn=rn*Xn*(1-Xn);
}

int iterations=lyapunovParams.iterations;
for (int n=0;n<iterations;n++){
 float prod_deriv=1;
 for (int i = 0; i < 3; i++) {
  float rn = i==2 ? p.x : i==1? p.y : p.z;
  Xn=rn*Xn*(1-Xn);
  prod_deriv *= rn*(1-2*Xn);
 }
 float q=log(fabs(prod_deriv));
 sum+=q;
}

return sum/(iterations*3);
}

This is the Kernel code for computing the positions of the evaluating points:

__global__ void kernelGeneration(float3 *ptr, int dim)
{
int xx=threadIdx.x + blockIdx.x*blockDim.x;
int yy=threadIdx.y + blockIdx.y*blockDim.y;
int offset=xx + yy*blockDim.x*gridDim.x;

int res=dim*dim;
int z=truncf(offset/res);
res=offset-z*res;
int y=truncf(res/dim);
int x=res-y*dim;

float fx = x/(float)dim;
float fy = y/(float)dim;
float fz = z/(float)dim;

ptr[offset].x=fx;
ptr[offset].y=fy;
ptr[offset].z=fz;
}

And finaly this is the Kernel code for computing the colors of those points:

__global__ void kernelComputation(float3 *ptrPos, float4 *ptrCol, int dim)
{
int xx=threadIdx.x + blockIdx.x*blockDim.x;
int yy=threadIdx.y + blockIdx.y*blockDim.y;
int offset=xx + yy*blockDim.x*gridDim.x;

float3 pos=ptrPos[offset];

pos.x*=lyapunovParams.maxX-lyapunovParams.minX;
pos.y*=lyapunovParams.maxY-lyapunovParams.minY;
pos.z*=lyapunovParams.maxZ-lyapunovParams.minZ;
pos.x+=lyapunovParams.minX;
pos.y+=lyapunovParams.minY;
pos.z+=lyapunovParams.minZ;

float lambda=evalLyap(pos);
float alpha=1;

if (lambda<-lyapunovParams.cut_off) {
 lambda+=lyapunovParams.cut_off;
 float l=-lambda*lyapunovParams.factor;
 if (l<0)l=0;
 if (l>1)l=1;
 ptrCol[offset].x=l;
 ptrCol[offset].y=l;
 ptrCol[offset].z=l;
 ptrCol[offset].w=l;
} else {
 ptrCol[offset].x=0;
 ptrCol[offset].y=0;
 ptrCol[offset].z=0;
 ptrCol[offset].w=0;
 ptrPos[offset].x=0;
 ptrPos[offset].y=0;
 ptrPos[offset].z=0;
}
}

Hoy voy a comentar un pequeño programa que no está totalmente relacionado con el tema del raytrace pero sí con el 3D.

Desde hace tiempo he intentado representar una versión tridimensional del fractal de Markus-Lyapunov. La versión 3D extiende el cómputo del exponente una dimensión más. Para realizarlo he utilizado CUDA y OpenGL.

Gracias a CUDA se puede realizar computación utilizando la GPU, lo cual, si se tiene una buena tarjeta de vídeo, puede hacer que la velocidad de cálculo sea bastante rápida.

La variación del algoritmo es sencilla. La secuencias contienen una letra más, en este caso C, que representa el punto en la coordenada Z (o 'c'). Todo lo demás es similar al algoritmo anterior.

Para realizar los renders utilizo dos VBOs: uno para almacenar posiciones y otro almacenar para colores. Las posiciones están repartidas en un volumen de 256x256x256 elementos. Para computar los colores de cada voxel utilizo dos kernels. Uno realiza el computo de las posiciones y las almacena en un VBO. El otro utiliza esas posiciones para realizar el computo de color, que es el exponente de lyapunov en ese punto. Para los puntos en los que el exponente es mayor que cero (o si se desea mayor que un valor de corte dado) el punto se lleva al origen para que no moleste. Para los puntos en que el exponente es menor que cero (o que el valor de cutoff) se pinta con un color proporcional al módulo de ese exponente.

Región espacial (0,4)-(0,4)-(0,4)

Región espacial (1,4)-(1,4)-(1,4)

Ajuste del valor de corte para quitar algo de información de estabilidad.

Reajuste del valor de corte para que esté más cercano a cero.

El código para computar el exponente de Lyapunov en 3D es este:

__device__ float evalLyap(float3 p)
{
float sum=0;
float Xn=0.5f;


for (int i = 0; i < 3; i++) {
 float rn = i==2 ? p.x : i==1? p.y : p.z;
 Xn=rn*Xn*(1-Xn);
}

int iterations=lyapunovParams.iterations;
for (int n=0;n<iterations;n++){
 float prod_deriv=1;
 for (int i = 0; i < 3; i++) {
  float rn = i==2 ? p.x : i==1? p.y : p.z;
  Xn=rn*Xn*(1-Xn);
  prod_deriv *= rn*(1-2*Xn);
 }
 float q=log(fabs(prod_deriv));
 sum+=q;
}

return sum/(iterations*3);
}

Este es el código de Kernel para computar las posiciones de los puntos de evaluación:

__global__ void kernelGeneration(float3 *ptr, int dim)
{
int xx=threadIdx.x + blockIdx.x*blockDim.x;
int yy=threadIdx.y + blockIdx.y*blockDim.y;
int offset=xx + yy*blockDim.x*gridDim.x;

int res=dim*dim;
int z=truncf(offset/res);
res=offset-z*res;
int y=truncf(res/dim);
int x=res-y*dim;

float fx = x/(float)dim;
float fy = y/(float)dim;
float fz = z/(float)dim;

ptr[offset].x=fx;
ptr[offset].y=fy;
ptr[offset].z=fz;
}

Y finalmente este es el código de kernel para computar los colores de dichos puntos:

__global__ void kernelComputation(float3 *ptrPos, float4 *ptrCol, int dim)
{
int xx=threadIdx.x + blockIdx.x*blockDim.x;
int yy=threadIdx.y + blockIdx.y*blockDim.y;
int offset=xx + yy*blockDim.x*gridDim.x;

float3 pos=ptrPos[offset];

pos.x*=lyapunovParams.maxX-lyapunovParams.minX;
pos.y*=lyapunovParams.maxY-lyapunovParams.minY;
pos.z*=lyapunovParams.maxZ-lyapunovParams.minZ;
pos.x+=lyapunovParams.minX;
pos.y+=lyapunovParams.minY;
pos.z+=lyapunovParams.minZ;

float lambda=evalLyap(pos);
float alpha=1;

if (lambda<-lyapunovParams.cut_off) {
 lambda+=lyapunovParams.cut_off;
 float l=-lambda*lyapunovParams.factor;
 if (l<0)l=0;
 if (l>1)l=1;
 ptrCol[offset].x=l;
 ptrCol[offset].y=l;
 ptrCol[offset].z=l;
 ptrCol[offset].w=l;
} else {
 ptrCol[offset].x=0;
 ptrCol[offset].y=0;
 ptrCol[offset].z=0;
 ptrCol[offset].w=0;
 ptrPos[offset].x=0;
 ptrPos[offset].y=0;
 ptrPos[offset].z=0;
}
}