julia set calculation with real numbers

Question

I am trying to make a julia set image and programming this in C . I have tried researching an algorithm to create an image of the Julia set but am getting confused by most of the examples online (as most examples seem to be copied but with little explanation of what happens step-by-step).

I have re-written my code several times to fit with algorithms I have found online but with little success. Currently, I am using an iteration function to define a colour of each pixel in the following:

// Loop through each pixel to work out new image set
int x, y;
for(x = 0; x < SIZE; x++){
    for(y = 0; y < SIZE; y++){

        int i = iterate(x, y, maxIterations, c);
        image[x][y] = i % 256 + 255 + 255 * (i < maxIterations);

    }
}

/* Iterate function */
int iterate(int x, int y, int maxI, double c[]){
    double z, zx, zy, oldRe, oldIm;   //real and imaginary parts of new and old
    double xmin = -1.0, xmax = 1.0, ymin = -1.0, ymax = 1.0;
    int k; // number of times iterated

    //calculate the initial real and imaginary part of z
    // z0 = (x + yi)^2 + c = (0 + 0i) + c = c
    zx = 1.5*(x - SIZE/2)/(0.5*SIZE);
    zy = 1.0*(y - SIZE/2)/(0.5*SIZE);

    //start the iteration process
    for(k = 1; k < maxI; k++){
      //remember value of previous iteration
      oldRe = zx;
      oldIm = zy;
      z = zx*zx - zy*zy + c[0];

      //the actual iteration, the real and imaginary part are calculated
      zx = oldRe * oldRe - oldIm * c[1] + c[0];
      zy = 2 * oldRe * oldIm + c[1];
      zy = 2.0*zx*zy + c[1];
      zx = z;

      //if the point is outside the circle with radius 2: stop
      if((zx * zx + zy * zy) > 4) break;
}

/*
while(((zx * zx + zy * zy) > 4) && (k < maxI)) {
    //remember value of previous iteration
    z = zx*zx - zy*zy + c[0];
    zy = 2*x*y-c[1];
    zx = z;
    //if the point is outside the circle with radius 2: stop
}
return k;

} */ The Commented out while loop is a second algorithm I found online but is also incorrect. I am not sure what I am missing from the equation. Any guesses?

(I am not looking for a solution, just some pseudocode that can convert the real, imaginary complex number to a position in my bounding box, in this case: -1,1).

complex number type is built in type in c , so one can use z = z*z + c — Adam, Nov 12 '17 at 19:10

score 0 · Answer 1 · answered Oct 31 '17 at 04:42

Here is the retranscription in readable code of a school project I just got validated

YOUR_CR and YOUR_CI are the two possible variables for Julia, LOWER_LIMIT and UPPER_LIMIT (X and Y) corresponds to the window you want to draw in. The code should scale correctly the output to the window and zoom level you want to draw the fractal in (the limits scale is in mathematical limits of the fractal, not in pixels).

# define X 0
# define Y 1

# define ZR 0
# define ZI 1
# define CR 2
# define CI 3

static int          iterations(long double coords[])
{
    long double     pos[4];
    int             i;
    double          tmp;

    pos[CR] = YOUR_CR;
    pos[CI] = YOUR_CI;
    pos[ZR] = coords[X];
    pos[ZI] = coords[Y];
    i = 0;
    while (i < MAX_ITERATIONS && (PYTHAGORE(pos[ZR], pos[ZI]) < 4))
    {
        tmp = pos[ZR];
        pos[ZR] = pos[ZR] * pos[ZR] - pos[ZI] * pos[ZI] + pos[CR];
        pos[ZI] = 2 * pos[ZI] * tmp + pos[CI];
        i++;
    }
    return (i);
}

# define TRANSLATE_X_PIXEL(a, b) ((int)(((a - LOWER_LIMIT_X) / (UPPER_LIMIT_X - LOWER_LIMIT_X)) * WIDTH))
# define TRANSLATE_Y_PIXEL(a, b) ((int)(((a - LOWER_LIMIT_Y) / (UPPER_LIMIT_Y - LOWER_LIMIT_Y)) * HEIGHT))

void                julia()
{
    long double     coords[2];
    int             itera;

    coords[X] = LOWER_LIMIT_X;
    while ((coords[X] += ((UPPER_LIMIT_X - LOWER_LIMIT_X) /
                    (double)(WIDTH + 1))) <= UPPER_LIMIT_X)
    {
        coords[Y] = LOWER_LIMIT_Y;
        while ((coords[Y] += ((UPPER_LIMIT_Y - LOWER_LIMIT_Y) /
                        (double)(HEIGHT + 1))) <= UPPER_LIMIT_Y)
        {
                itera = iterations(coords);
                set_pixel(TRANSLATE_X_PIXEL(coords[Y]), TRANSLATE_Y_PIXEL(coords[Y]), itera * DESIRED_COLOR);
        }
    }
}

score 0 · Answer 2 · answered Oct 31 '17 at 07:40

The out-commented while code is wrong as it uses x,y in the imaginary part computation instead of zx,zy. Your code looks as if you took one approach and pasted another in the middle of it.

What you want the classical quadratic Julia iteration, it is the realization in real components of

znext = z*z + c = (zx+*zy)*(zx+i*zy)+(cx+*cy)
      = (zx*zx-zy*zy+cx) + i*(2*zx*zy+cy)

You can not do in-place replacements since both zx and zy are used in both formulas and need to keep their old values there. The correct code should be

for(k = 1; k < maxI; k++){
    // compute the new real part, store temporarily in helper variable
    double z = zx*zx - zy*zy + c[0];
    // compute the imaginary part, since no further operations, do it in place
    zy = 2.0*zx*zy + c[1];
    //store the new real part in the real part variable
    zx = z;
    //if the point is outside the circle with radius 2: stop
    if((zx * zx + zy * zy) > 4) break;
}

You would need less multplications if you store the squares,

for(k = 1; k < maxI; k++){
    double sqzx = zx*zx, sqzy = zy*zy;
    // compute the imaginary part first, since the real part only uses the squares, not zy
    zy = 2.0*zx*zy + c[1];
    // compute the new real part using the pre-computed squares
    zx = sqzx - sqzy + c[0];
    //if the point is outside the circle with radius 2: stop
    if( ( sqzx + sqzy ) > 4 ) break;
}

julia set calculation with real numbers

2 Answers2