Stripe Average Coloring

[asimes]

I found a thesis from a post here that explains a number of fractal coloring techniques by Jussi Harkonen. The thesis can be downloaded here: http://jussiharkonen.com/?page_id=65

He gives an explanation on how to do Stripe Average Coloring but to be honest I can't understand half of his equations (I never took Calculus or Complex Analysis). Does anyone know of another resource to learn how to do this? Google didn't seem too helpful

Incase it isn't clear what kind of coloring technique I mean, I mean something like this: http://www.fractalforums.com/2d-art/stripe-average-mandelbrot-2/

[Syntopia]

Nice find. I tried implementing the Stripe Average scheme and it is a great coloring function.

The basics steps are these:
1) Calculate the orbit for a given point in the complex plane. This is series of complex numbers {z0,z1,z2,...}
2) Apply some function to this series to get a new series of real numbers: { f(z0), f(z1), f(z2), ...}. His 'stripe' function is basically f(z)=sin(arg(z)) plus some scalings.
3) Take the average of this series of real numbers.

This gives you a real number, which you have to map somehow to a palette. If you do this, you get images like the first attached one.

Now things get a little tricky, because of some discontinuities whenever there is a jump in the integral number of iterations before the orbit escapes. But it is possible to construct a 'smooth' fractional iteration count, and he uses the fractional part of this to interpolate between the color value calculated from the original series and the color value obtained from the same series with the last element excluded.

The fractional part of the smooth iteration count is:

Code:

float frac =-1.0+log2(.5*log(EscapeRadius*EscapeRadius))-log2(.5*log(lastZ2));

where lastZ2 is the length of the last point in the orbit, before it escapes.

The second image shows how this removes the iteration discontinuities. There are still some discontunuities, but these can be removed by skipping a fixed number of items from the sequence.

I've attached a script with some Fragmentarium code that might help.

[asimes]

Ok, I think I finally made sense of how to get average color working. I attached an image of the Mandelbrot that I think has the correct stripes below.

Unfortunately, I don't understand how to smooth it out. I have a few questions about it:

1. What do I do with that scary looking float, frac?
2. How do I store lastZ2? Is it just the last iteration of sqrt(zReal*zReal+zImaginary*zImaginary) given a point?
3. In my code I decided to get a real number sequence by using sin(zReal+zImaginary) in a while loop. Is this the same as f(z) = sin(arg(z))?

Thanks for the help

[asimes]

I just noticed that my stripes are arcs and yours are lines perpendicular to the iteration bands. I think I did something wrong but I don't know what. I wasn't sure on how to write out the orbit in code so I just made this (in pseudo-code):

for (all of the points on the complex plane) {
  float zr = realCoordinate;
  float zi = imaginaryCoordinate;
  float orbitCount = 0;
  int n = 0;
  while (n < maxIterations) {
    float zrr = zr*zr;
    float zii = zi*zi;
    float twori = 2*zr*zi;
    orbitCount += sin(zr+zi);
    zr = zrr-zii+x;
    zi = twori+y;
    if (zrr+zii > 16.0) break;
    n++;
  }
  if (n == maxIterations) make the point black
  else {
    orbitCount /= n;
    float orbitColor = (orbitCount+1)*127.5;
    make the point orbitColor
  }
}

[Syntopia]

Ok, I think I finally made sense of how to get average color working. I attached an image of the Mandelbrot that I think has the correct stripes below.

Unfortunately, I don't understand how to smooth it out. I have a few questions about it:

1. What do I do with that scary looking float, frac?
2. How do I store lastZ2? Is it just the last iteration of sqrt(zReal*zReal+zImaginary*zImaginary) given a point?
3. In my code I decided to get a real number sequence by using sin(zReal+zImaginary) in a while loop. Is this the same as f(z) = sin(arg(z))?

Thanks for the help

<Quoted Image Removed>

1) You need it to interpolate between two values, one obtained from the full series, and one obtained from the average of all elements except the last one:
float mix = frac*avg+(1.0-frac)*prevAvg;

2) Yes, the last orbit value with a length less then the escape radius. Notice, lastZ2 is the length squared, so don't take the square root.
   
3) No, you need the argument, which is the angle of the complex vector to the real axis. I used this function:
f(z) = 0.5+0.5*sin(StripeDensity*atan(z.y,z.x))

[DarkBeam]

It's plenty of working coloring algorithms in the public ultra fractal database (thousands) , no need to break your head to find new

[asimes]

Well, I got the stripes working for real this time (my Julia Set version looks like the Julia Set from the thesis paper examples), I'll post a picture of it below. I still don't understand something from your new explanation I think. I tried adding the new information you mentioned but the result still has iteration bands. Here is the code if you have the patience to look through it, although I feel like the problem comes from not know what to do with the float mix:

Code:

float xmin = -2.5; 
float ymin = -2.0; 
float wh = 4;
float stripes = 5.0;
int maxIterations = 1000;

void setup() {
  size(800, 800, P2D);
}

void draw() {
  loadPixels();
  float xmax = xmin+wh;
  float ymax = ymin+wh;
  float dx = (xmax-xmin)/width;
  float dy = (ymax-ymin)/height;
  float x = xmin;
  for (int i = 0; i < width; i++) {
    float y = ymin;
    for (int j = 0;  j < height; j++) {
      float zr = x;
      float zi = y;
      float lastzr = x;
      float lastzi = y;
      float orbitCount = 0;
      int n = 0;
      while (n < maxIterations) {
        float zrr = zr*zr;
        float zii = zi*zi;
        float twori = 2*zr*zi;
        zr = zrr-zii+x;
        zi = twori+y;
        if (zrr+zii > 16) break;
        orbitCount += 0.5+0.5*sin(stripes*atan2(zi, zr));
        lastzr = zr;
        lastzi = zi;
        n++;
      }
      if (n == maxIterations) pixels[i+j*width] = 0;
      else {
        float lastOrbit = 0.5+0.5*sin(stripes*atan2(lastzi, lastzr));
        float smallCount = orbitCount-lastOrbit;
        orbitCount /= n;
        smallCount /= n-1;
        float frac = -1+log(0.5*log(16))-log(0.5*log(lastzr*lastzr+lastzi+lastzi));
        float mix = frac*orbitCount+(1-frac)*smallCount;
        float orbitColor = mix*255;
        pixels[i+j*width] = color(orbitColor);
      }
      y += dy;
    }
    x += dx;
  }
  updatePixels();
  noLoop();
}

Here is the image of the new stripes for the non-smoothed average color stripes (not my attempt at smoothing which looks almost the same really at the moment):

[asimes]

@DarkBeam, I'm not familiar with it, but I'm trying to learn how to code them myself to make sense of it (obviously it is still over my head).

[Syntopia]

Hmm, I think the problem might be in your 'frac' line.

Code:

  float frac = -1+log(0.5*log(16))-log(0.5*log(lastzr*lastzr+lastzi+lastzi));      

Try using log-2 instead of log-10 for the outer logarihms. If you haven't got them, divide by log(2):

Code:

  float frac = -1+(log(0.5*log(16))-log(0.5*log(lastzr*lastzr+lastzi+lastzi)))/log(2);

Alnd try raising the escape radius to something much larger than 4.

The difficult part is using the right orbit point in the formulas: actually, I think I have an error my original code - the length inside the logarithm in the line above, should probably be evaluated at the second last element, instead of the last (and the last element should be removed from both series). You don't have to fix this, though. Just try changing the coefficient in front of log(16) from 0.5 to 1.0 should have the same effect.

Attached a Mandelbrot example.

[asimes]

This is what I've got for log(): http://processing.org/reference/log_.html

I tried playing with the numbers and 'frac' but I can't seem to make the iterations bands disappear so I think I'm not getting something fundamental. It definitely looks nicer but none of the changes I make affect the iteration bands, just the contrast of the stripes and their width / density. How are you making your images? Would you be willing to post your code? I'm starting to feel like I'm asking too much to keep posting but at least it is close (ignoring the iteration bands):

Code:

float xmin = -2.5; 
float ymin = -2.0; 
float wh = 4;
float stripes = 5.0;
int maxIterations = 1000;

void setup() {
  size(800, 800, P2D);
}

void draw() {
  loadPixels();
  float xmax = xmin+wh;
  float ymax = ymin+wh;
  float dx = (xmax-xmin)/width;
  float dy = (ymax-ymin)/height;
  float x = xmin;
  for (int i = 0; i < width; i++) {
    float y = ymin;
    for (int j = 0;  j < height; j++) {
      float zr = x;
      float zi = y;
      float lastzr = x;
      float lastzi = y;
      float orbitCount = 0;
      int n = 0;
      while (n < maxIterations) {
        float zrr = zr*zr;
        float zii = zi*zi;
        float twori = 2*zr*zi;
        zr = zrr-zii+x;
        zi = twori+y;
        if (zrr+zii > 10000) break;
        orbitCount += 0.5+0.5*sin(stripes*atan2(zi, zr));
        lastzr = zr;
        lastzi = zi;
        n++;
      }
      if (n == maxIterations) pixels[i+j*width] = 0;
      else {
        float lastOrbit = 0.5+0.5*sin(stripes*atan2(lastzi, lastzr));
        float smallCount = orbitCount-lastOrbit;
        orbitCount /= n;
        smallCount /= n-1;
        float frac = -1+log(0.5*log(10000))/log(2)-log(0.5*log(lastzr*lastzr+lastzi*lastzi))/log(2);
        float mix = frac*orbitCount+(1-frac)*smallCount;
        float orbitColor = mix*255;
        pixels[i+j*width] = color(orbitColor);
      }
      y += dy;
    }
    x += dx;
  }
  updatePixels();
  noLoop();
}

[Syntopia]

The code is posted as an attachment to my first post - it is the ".frag" file. It is GLSL, but should be readable. Ignore the map-color functions and presets - the code is in the getColor2D-function.

I think your log function should be okay, even if it is base-e - as long as you divide by log_e(2), you should get a log_2() functions.

[Syntopia]

I just tested your code (I realized I had processing installed). If you change to

Code:

 float frac = -1+log(2.0*log(10000))/log(2)-log(0.5*log(lastzr*lastzr+lastzi*lastzi))/log(2);      

it should work. (It is difficult to explain, but the square-length of your 'lastz' variable may be larger than the squared escape distance, the way you setup your break condition - and we really needed the variable two iterations earlier. However, for large z, we have quadratic growth, so we can just square the escape radius twice - corresponding to multiplying by 4 outside the logarithm).

[asimes]

That worked! Thank you

To be honest, I don't know why changing 0.5 to 2.0 made it work, but it looks great now. I was looking through your code (I didn't realize you had posted it until you mentioned it in your last response) and I liked how you did the orbit adding so I'm going to change mine.

[Syntopia]

The

That worked! Thank you

To be honest, I don't know why changing 0.5 to 2.0 made it work, but it looks great now. I was looking through your code (I didn't realize you had posted it until you mentioned it in your last response) and I liked how you did the orbit adding so I'm going to change mine.

The expression got unnecessarily ugly, because I choose the wrong escape point. The correct expression should be:

Code:

float frac =1.0+(log2(log(EscapeRadiusSquared)/log(z2)));

where z2 is the square length of the first orbit point outside the escape radius. For some reason I always mange to offset these orbit :-)

A nicer version of the script is added.

Btw, In the thesis he uses spline interpolation instead of linear interpolation, which should give even better results. But I haven't tried this yet.

[asimes]

Well, I know enough about manipulating the log() at this point at least to convert your new one to one Processing can use. The Processing usable one is this now (because there is no log2):

float frac = 1+(log(log(10000)/log(lastLength))/log(2));

I don't understand enough about what is going on conceptually to be able to attempt a spline interpolation because I don't know what frac's is doing. I understand what the rest of the code does but frac is confusing to me (this is also the first time I make serious use of log in a sketch).