Helix fractals using raymarching

The goal of this article is to give you some ideas on how to create fun helix fractals in GLSL. Recently I made some gif animations using this algorithms:

Untangle fractal 1.02

Untangle fractal 1.50

Untangle fractal 2.00

These pieces may seem complex at the first sight, but they all share a really simple algorithm, which you may probably guess what it is if you pick up a pen and paper and spent a few minutes thinking about it. To start we need a simple raymarching code which can render a simple sphere for you. It'll be good if your code has soft-shadow and ambient occlusion already implemented.

The essence of the raymarching algorithm is to have a "Shortest Distance Fields" SDF(pos) which in every point of space gives you the shortest distance to the surface that light rays may hit. We can render all sort of complex objects just by finding a way for approximating the SDF, or at least a function that we know is always <= to the real SDF.

First Draw a cylinder, The SDF is pretty easy in this case:

float SDCyl(vec3 pos,float R){ return length(pos.xz) - R; }

figure 1 : Raw cylinder


Now we must pass the pos vector from some mapping function before passing it to the cylinder SDF. What this mapping should look a like? 

Our goal is to define a function that takes a cylinder and remaps it to a helix. We call it "twist"

figure 2


First step: copy the cylinder
    
figure 3

as you can see from the above picture we would like to cut the space into 4 parts (or more generally N pieces). and in each part of the space we define a new center of position and we rotate the x and z vectors into x' and  z'. For now consider N=4. Coding this would be easy:

vec3 twist(vec3 pos, float t , float R ){
    // "t" is the angle offset defined in figure 3
    //cylindrical coordinate
    float theta = atan(pos.z,pos.x);
    float r = length(pos.xz);

    //Just bringing theta to the desired interval and decreasing the offset value "t" from it
    theta = mod(theta - t , 2.*PI);

    //Using round(float) function to find the theta of the nearest circle
    float nt = round(theta/(2.*PI/N))*2.*PI/N;
    nt += t; //don't remember the offset must be added at the end

    //now we find the center position of the nearest circle
    vec3 np = vec3(R*cos(nt) , pos.y , R*sin(nt));

    //finding x',y',z' by rotating x,y,z
    vec3 yp = vec3(0. , 1. , 0.);
    vec3 xp = vec3(cos(nt),0.,sin(nt));
    vec3 zp = cross(xp,yp);

    //the pos we gave to the twist function now must be measured in the new coordinate system
    vec3 ret = vec3(dot(pos-np,xp) , pos.y, dot(pos-np,zp));    

    //return the mapped position
    return ret; 
}


Now let's check what this SDF would be like: 



figure 4: SDCyl(twist(pos,0.,1.),0.5);

As you can see the code copies the cylinder 4 times and puts them on a circle, in order to do that, for every point in space, it calculates the position in cylindrical coordinate, r is distance from the y axis, and theta is the azimuthal angle. Then it check which of the circles in figure 3 is the nearest, but before that we use the "t" input as an offset.  We call the nearest circle's theta "nt" can be found using round() function.

Now what if we rotate the offset value t as we go up on the cylinder? This means setting:

t += pos.y/R;

at the beginning of the twist function. Then we suddenly get a helix like this:

figure 5: A nice helix, but something is wrong.

Do you see that too? There is something off about this helix! It's like that the helix is extra thin vertically! You may guess what is the problem, I'll give you a hint, let's color the helix based on the y' value of each point and see how it looks:

figure 6: The coloring hints us the problem

Now as you may guess the problem is  y'=y ! this must not be the case for a nice round helix! The x' y' z' coordinates must now change in a tricky way:

   vec3 yp = vec3(-sin(nt) , 1. , cos(nt))/sqrt(2.); //y' must be rotated around x'
   vec3 xp = vec3(cos(nt),0.,sin(nt));
   vec3 zp = cross(xp,yp);
   ret = vec3(dot(pos-np,xp) , pos.y*sqrt(2.)  + dot(pos-np,yp), dot(pos-np,zp));   //Think about the piece of code added in this line.

Then the result would be:

figure 7: Smoooooth!

We are now in a very good place. Let's do something fun. Before calling the SDF function pass the position through this function:

    float ALPHA = 0.30; // The scaling ratio for each iteration
    float r_small = 1.5; // The radius of the cylinder after all the copying and twisting
    float R = 1.5; //Radius of the copying circle
  
    for(int i=0; i<4;i++){
        pos  = twist(pos, PI/8. , R,0.) ;
        r_small *= ALPHA;
        R *= ALPHA;
    }

Now let's call the SDF and see what does it draw for us:

figure 8: Calling twist 4 times, makes a helix in helix fractal

This shape is one of my most favorite fractals! You can also find it in nature, in the structure of DNAs! Check out this awesome video.

Let's test some other cool effect. Let's make the radius R in twist function larger in some places. Add this in the beginning of the twist function:

R = R*(1. + 10.*(0.5+0.5*cos(clamp(pos.y/R/60.,-1.,1.)*PI)));

Then zoom out and watch the result!

figure 9: This is how I made "Untangled fractal 2.0"

Play around with what we got until here. I also made this: (using another function to inflate R and also turning off the rotation of offset "t"):


figure 10: another one!


Another thing we could do: In the twist function, do different mappings for pos.y>0 and pos.y<0.
And try to make it like this. Also offset pos.y in each iteration by a value of  -H*ALPHA^i. 
figure 11: a new trick

Then if you use twist iteratively you'll get this:


figure 12: The result


Let me know what other ideas you've got. Please notify me if you you found a mistake in my code and algorithm or you know any different approaches. I would be more than happy to hear from you guys! 




   

Comments

Post a Comment