Sunday, April 12, 2015

Mandelbrot Set Rendering

Benoit Mandelbrot was a Polish-born mathematician. His most well known work is in the field of fractals. He came up with a function describing a set of complex numbers called the Mandelbrot Set.
The Mandlebrot set describes a set of numbers "c" such that the sequence (c1, c2 .... cx) where cn+1 = (cn*cn)+cn does not approach infinity.
So you're looking at c, c²+c, (c²+c)²+c, ((c²+c)²+c)²+c...

A complex number has two parts. It's a sum of real an imaginary components. Let's say we're looking at A + Bi.
A is the real component.
Bi is the imaginary component.
'i' is a symbol representing the √-1 (square root of -1)
i.e., (Bi)² would be -B².

Mandelbrot's equation for describing the model is z'=z²+c.
Let z be our A + Bi from above.
The next element of the sequence is 
If z=(A+Bi)² then z=A²-B²+2*A*Bi


If plotting, we can us X, Y in place of A, and B.
so, z=X²-Y²+2*X*Y(i)
We have real and imaginary components here.
So Xnew=X²-Y², and Ynew=2*X*Y
the "c" is the starting point for Mandelbrot's equation.
Z=Z²+c, so we need to add X0 and Y0 to the above:

Xnew=X²-Y²+X0
Ynew=2*X*Y+Y0
And from here, we keep iterating on this till we hit set limits to number of iterations, or if X/Y are sufficiently large.

So, for a point X, Y,
Repeat the following until x0²+y0²>4 or until set max iterations

Until (max iterations, or x0²+y0²>4)
  x0=0 and y0=0
  x1=x0*x0-y0*y0+X
  y1=2*x0*y0+Y
  x0=x1
  y0=y1
Color point X,Y based on iterations.




Point X,Y is then associated with the number of iterations to break out of that loop.

Using this we can construct a basic Mandelbrot generator. I did this using Javascript, since it's easy to use any browser to display a canvas element. This works in Firefox and Chrome (has worked for a while) and I've tested it in IE11 and it works there too.




<!DOCTYPE HTML>
<html>
<script language="javascript">
  function rgbToHex(r, g, b) { //thanks stackoverflow
    return "#" + ((1 << 24) + (r << 16) + (g << 8) + b).toString(16).slice(1);
  }
  function putpixel(x,y,r,g,b,ctx){
    ctx.fillStyle = rgbToHex(r,g,b);
    ctx.fillRect(x,y,1,1);
  }
  function paint(cnvname){
    mcanvas = document.getElementById(cnvname);
    ctx = mcanvas.getContext('2d');
    ctx.clearRect(0, 0, 799, 599);
    for (x = 0; x < 799; x++) {
      for (y = 0; y < 599; y++) {
        tx=(x-400)/(200); // tx is "translated x" for graph centering"
        ty=(300-y)/(200); // ty is "translated y" for graph centering"
        xtmp=0;
        ytmp=0;
        xtmp2=0;
        iter=0
        maxiter=1000;
        for (; (iter<maxiter) && (xtmp*xtmp+ytmp*ytmp < 4);) {
          xtmp2=xtmp*xtmp-ytmp*ytmp+tx;
          ytmp=2*xtmp*ytmp+ty;
          xtmp=xtmp2;
          iter++;
        }
        contrast=0.15
        color=Math.floor( (Math.pow(iter,contrast)*255)/Math.pow(maxiter,contrast));
        putpixel(x,y,0,255-color,255-color,ctx);
      }
    }
  }
</script>
<body onload='paint("mcanvas")'>
  <div>
    <canvas id="mcanvas" width="800" height="600"></canvas>
  </div>
</body>
</html>


This is a simple bit of javascript that iterates over the pixels of the canvas created in the body. It applies the function described earlier to determine the max number of iterations to the set.



Pretty right? Of course there's much more to mandelbrot than just the bulb.
We can zoom in on a portion of this by applying a zoom and x/y offset to the above code.

We can add a zoom and offset feature to explore the set.


  for (x = 0; x < 799; x++) {
    for (y = 0; y < 599; y++) {
      zoom=20; xoff=0.6; yoff=0.6;
      tx=-xoff+(x-400)/(zoom*200); // tx is "translated x" for graph centering"
      ty=yoff+(300-y)/(zoom*200); // ty is "translated y" for graph centering"

The above change will produce the following:

This is zoomed in 20x onto the small bulb in the top left portion of the largest bulb. Note the pattern replicates, and there are several tiny bulbs that look a lot like the main one (with slight variations).

Now the fun part - we can color based on x, y or iterations. Apply a sinusoidal function to iterations for red, blue and green, and you can have a never ending palette of color.

e.g.
From the contrast and color line:

//contrast=0.25
//color=Math.floor( (Math.pow(iter,contrast)*255)/Math.pow(maxiter,contrast));
cycle=(maxiter)/(3.14159265359*100);
cr=Math.floor(128+127*Math.sin(0.00+iter/(cycle*1.5000)));
cg=Math.floor(128+127*Math.sin(1.57+iter/(cycle*1.6000)));
cb=Math.floor(128+127*Math.sin(3.14+iter/(cycle*1.5666)));
putpixel(x,y,cr,cg,cb,ctx);

So, what are we saying? A cycle is 2*pi. We're dividing the iterations by 100*pi, so we're looking at at 50 cycles. Red, green and blue values are sinusoidal cycles with different start points - red starts at 128+127*sin(0) (so 128, and increasing). Green is 128+127*sin(1.57) (so 255, 1.57 being half pi, sin(pi/2) is 1) - i.e. max green and heading down to 0. Blue is starting at 128+sin(pi), which means middle intensity and getting darker.
Each sinusoidal cycle has it's own multiplier so each cycle has it's own rate, giving many more colors ans each cycle creates a new palette.

Here's the result:


We're looking at integers since we're only counting max iterations. We can smoothen this by looking at how far past our boundary condition each iteration went.

My script was based on methods described in wikipedia.


log2=Math.log(2)
if (iter<maxiter){ //smoother transitions
  zn=Math.sqrt(xtmp*xtmp+ytmp*ytmp);
  zl=Math.log(Math.log(zn)/log2)/log2;
  iter+=1-zl;
} 






Now consider the interior of the bulb. It's what happens when the iterations have been maxed. We can color the interior by using the same palette. I haven't quite figured out the how to map distance, but we can use the existing xtmp and ytmp variables. Since iterating more doesn't get us greater than the boundary condition, we can assume xtmp and ytmp are small. 

What I tried was several permutations of setting iter to 1/xtmp and 1/ytmp. I settled on using iter=Math.log( (1/(xtmp*xtmp+ytmp*ytmp)) )/(log2);
The idea is that if my values aren't going to break out of the boundary, then they're converging. So I'm using the inverse of the values to dictate what they will be against the existing palette. Not all sub-bulbs show this though.


Here's the overall code for this basic mandelbrot renderer...



<!DOCTYPE HTML>
<html>
<script language="javascript">
  function rgbToHex(r, g, b) { //thanks stackoverflow
    return "#" + ((1 << 24) + (r << 16) + (g << 8) + b).toString(16).slice(1);
  }
  function putpixel(x,y,r,g,b,ctx){
    ctx.fillStyle = rgbToHex(r,g,b);
    ctx.fillRect(x,y,1,1);
  }
  function paint(cnvname){
    mcanvas = document.getElementById(cnvname);
    ctx = mcanvas.getContext('2d');
    ctx.clearRect(0, 0, 799, 599);
    log2=Math.log(2)
    for (x = 0; x < 799; x++) {
      for (y = 0; y < 599; y++) {
        zoom=20; xoff=0.6; yoff=0.6;
        tx=-xoff+(x-400)/(zoom*200); // tx is "translated x" for graph centering"
        ty=yoff+(300-y)/(zoom*200); // ty is "translated y" for graph centering"
        xtmp=0;
        ytmp=0;
        xtmp2=0;
        iter=0
        maxiter=1000; // increase with larger zoom values
        totes=0;
        for (; (iter<maxiter) && (xtmp*xtmp+ytmp*ytmp < 4);) {
          xtmp2=xtmp*xtmp-ytmp*ytmp+tx;
          ytmp=2*xtmp*ytmp+ty;
          xtmp=xtmp2;
          iter++;
        }
        if (iter<maxiter){ //smoother transitions
          zn=Math.sqrt(xtmp*xtmp+ytmp*ytmp);
          zl=Math.log(Math.log(zn)/log2)/log2;
          iter+=1-zl;
        } else { // for internal color of Mandelbulb
          iter=Math.log( (1/(xtmp*xtmp+ytmp*ytmp)) )/(log2);
        }
        cycle=(maxiter)/(3.14159265359*100); //50 cycles from 0-maxiter
        cr=Math.floor(128+127*Math.sin(0.00+iter/(cycle*1.5000)));
        cg=Math.floor(128+127*Math.sin(1.57+iter/(cycle*1.6000)));
        cb=Math.floor(128+127*Math.sin(3.14+iter/(cycle*1.5666)));
        putpixel(x,y,cr,cg,cb,ctx);
      }
    }
  }
</script>
<body onload='paint("mcanvas")'>
  <div>
    <canvas id="mcanvas" width="800" height="600"></canvas>
  </div>
</body>
</html>


Manipulation of the color palette constants, zoom, offsets, and even the equation can produce some beautiul results.

For example,
Set the initial xtmp to tx and ytmp to ty.

change:
xtmp2=xtmp*xtmp-ytmp*ytmp+tx;
ytmp=2*xtmp*ytmp+ty;


to:
xtmp2=xtmp*xtmp-ytmp*ytmp-0.32;
ytmp=2*xtmp*ytmp+0.6575;


Now you're making a Julia plot.
This is another fractal you can zoom in continuously.

Let's look at zooming in on a portion of the Mandelbrot set.
zoom=1.00; xoff=0; yoff=0;


zoom=5; xoff=0.7348; yoff=0.1798;

zoom=25; xoff=0.8; yoff=0.2;



zoom=125; xoff=0.79; yoff=0.158;


zoom=625; xoff=0.79; yoff=0.162;

zoom=3125; xoff=0.7893; yoff=0.163;


At this point we need to increase the maxiter constant. I raised this to 10000.
zoom=9125; xoff=0.78937; yoff=0.16307;


 zoom=15625; xoff=0.78940; yoff=0.163055;


zoom=78225; xoff=0.78940; yoff=0.163055;
Note the zoom - we're at almost 80000x and as you can see there's a lot more structures here!