Julia set

  • 13 Replies
  • 400 Views

0 Members and 1 Guest are viewing this topic.

Offline Adam Majewski

  • *
  • Fractal Friar
  • *
  • Posts: 126
« on: December 29, 2018, 09:52:24 AM »
Hi,

How these images :
http://virtualmathmuseum.org/Fractal/julia/index.html
were created ?

What algorithm was used ?
I mean : julia 01  and next where parts of Julia set have different color

TIA

Adam

Offline FractalDave

  • *
  • Uploader
  • *
  • Posts: 171
    • Makin Magic Fractals
« Reply #1 on: December 31, 2018, 09:54:58 PM »
Colours: Section of genetic code based on binary for + or - square root on successive iterations ?
The meaning and purpose of life is to give life purpose and meaning.

Offline knighty

  • *
  • Fractal Feline
  • **
  • Posts: 185
« Reply #2 on: January 02, 2019, 12:54:02 PM »
I'm also curious about how it works.
It seems to be some kind of IIM with a choice of square root branch based on a "angle like" state.

Offline Adam Majewski

  • *
  • Fractal Friar
  • *
  • Posts: 126
« Reply #3 on: January 02, 2019, 08:41:35 PM »
Thx for the tips.

maybe this will explain smth :
https://arxiv.org/abs/hep-th/0501235
page 67, 75
?

Offline Softology

  • *
  • Fractal Fanatic
  • ***
  • Posts: 38
« Reply #4 on: January 02, 2019, 09:24:40 PM »
What algorithm was used ?
I mean : julia 01  and next where parts of Julia set have different color

Ask him.  Contact email is at the bottom of this page http://virtualmathmuseum.org/Fractal/index.html
Author home page http://vmm.math.uci.edu/

Offline knighty

  • *
  • Fractal Feline
  • **
  • Posts: 185
« Reply #5 on: January 03, 2019, 11:28:06 AM »
Failed attempt  :embarrass:

Still, it sometimes give correct result :


Anyway... The author gives a hint in the .pdf :
Quote
In  our  program  we  compute  preimages  starting  from the circle {z; |z|=Rc} around the wanted Julia set.
Under inverse images these curves converge from outside to the Julia set.  Such an approximation by curves allows  us  to  color  the  Julia  set  in  a  continuous  way and  thus  emphasize  that,  despite  its  wild  looks  it  is the image of a continuous curveat least for c ∈ M.

Offline Adam Majewski

  • *
  • Fractal Friar
  • *
  • Posts: 126
« Reply #6 on: January 03, 2019, 04:39:49 PM »
Ask him.  Contact email is at the bottom of this page http://virtualmathmuseum.org/Fractal/index.html
Author home page http://vmm.math.uci.edu/
I did it. No response

Offline Adam Majewski

  • *
  • Fractal Friar
  • *
  • Posts: 126
« Reply #7 on: January 03, 2019, 04:42:50 PM »
Failed attempt  :embarrass:

looks great, thx

So the algorithm is intended to show Jc as a continuously parametrized curve.

Steps
  • choose c
  • compute radius R = 1/2 + sqrt( 1/4 + |c|)
  • compute n points z = R * e^(i*t) from circle with center at origin and radius R. Color of point z is proportional to angle t   
  • for each point z make some inverse iteration  until it fall into the Julia set
  • color of the point z from Julia set is the same as color of point z from initial circle

I think that most important is to choose such backward iteration that point moves along the same  external ray (with external angle = t) toward Julia set.
It is described here as : Backward iteration (with proper chose from two possibilities)
https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/boettcher#backwards_iteration

angle can be taken as an rational number t = p/q to make iterations easier

Or you can use any other method for computing external ray  ( Newton, inverse Boettcher map)  to find it's landing point


So the method can look like this:


 

- it is a from https://gitlab.com/adammajewski/periodic-points-of-complex-quadratic-polynomial-using-newton-method

« Last Edit: January 06, 2019, 09:05:34 AM by Adam Majewski, Reason: new image »

Offline FractalDave

  • *
  • Uploader
  • *
  • Posts: 171
    • Makin Magic Fractals
« Reply #8 on: January 03, 2019, 08:25:10 PM »
Failed attempt  :embarrass:

Still, it sometimes give correct result :


Anyway... The author gives a hint in the .pdf :

OK, from an old wip formula for UF (use Basic colouring from standard.ucl set to real with linear transfer density the same as the iteration count in the parameters i.e. detail - I have no recollection of why I called the formula "Two", I'm guessing it's based on similar info to that the author used in the images this thread refers to):

Code: [Select]

MMF-Two {
global:
;  $define debug
;;  print((cos(135*#pi/180)+flip(sin(135*#pi/180)))/sqrt(2))
;;  print((cos(225*#pi/180)+flip(sin(225*#pi/180)))/sqrt(2))
  int i=0
  int j
  int k
  int l
  int m
  int n
  int s[@width+1,@height+1]
  complex v[2,round(2^@detail)];round(2^(24-@detail))]
  int angles = round(2^@angles)
  float q[282]
  float a
  complex c
  complex zz
  repeat
    a = i*#pi*2.0/angles
    zz = v[0,0] = 2.0*(cos(a)+flip(sin(a)))
    j = 0
    m = 0
    n = 1
    repeat
      k = 0
      l = round(2^j)
      repeat
        c = v[n,l] = sqrt(v[m,k]-@const) ;1.0-(0.5-flip(0.5))*v[m,k,i];0.5*(sqrt(1+4.0*v[m,k,i])-1.0);
        if @alliters || j==@detail-1
          s[floor(0.5*@width+real(c)*@width*#magn/4.0),floor(0.5*@height+imag(c)*@height*#magn/3.0)]=(l-round(2^j))*angles+i+1
        endif
        c = v[n,k] = -v[n,l] ;-1.0-(0.5+flip(0.5))*v[m,k,i];-0.5*(sqrt(1+4.0*v[m,k,i])+1.0);
        if @alliters || j==@detail-1
          s[floor(0.5*@width+real(c)*@width*#magn/4.0),floor(0.5*@height+imag(c)*@height*#magn/3.0)]=(l-round(2^j))*angles+i+1
        endif
        l = l + 1
      until (k=k+1)>=round(2^j)
      m = (m+1)%2
      n = (n+1)%2
    until (j=j+1)>=@detail
    ;print(i,"/",angles)
  until (i=i+1)>=angles
init:
  z = 0
  bool res = false
  float t
  if s[round(real(#screenpixel)),round(imag(#screenpixel))]>0
    res=true
    t = (s[round(real(#screenpixel)),round(imag(#screenpixel))]-1)
    z = (t%angles)/angles + flip((t-t%angles)/(angles*2^(@detail-1)))
  endif
loop:
bailout:
  !res
default:
  title = "Two"
  render=false
  complex param const
    default = (-0.12, 0.8)
    hint = "Julia seed"
  endparam
  bool param alliters
    default = false
  endparam
  int param detail
    default = 12
    max = 25
    min = 1
    hint = "max = 25 iterations"
  endparam
  int param angles
    caption = "#angles (2^)"
    default = 9
    min = 0
    max = 10
    hint = "max = 10 i.e. 2^10=1024"
  endparam
  int param width
    default = 640
    hint = "Should match image widrh"
  endparam
  int param height
    default = 480
    hint = "Should match image height"
  endparam
}


Offline FractalDave

  • *
  • Uploader
  • *
  • Posts: 171
    • Makin Magic Fractals
« Reply #9 on: January 03, 2019, 08:50:06 PM »
Or you can use any other method for computing external ray  ( Newton, inverse Boettcher map)  to find it's landing point

The brute force UF Wip formula in my last post was done before I knew about Boettcher, which means it's definitely pre 2008 ;)
In 2008 I used ideas related to the method for image mapping to external angles/smooth iteration:

https://www.deviantart.com/makinmagic/art/Fractals-are-hard-to-swallow-102574524

Offline Adam Majewski

  • *
  • Fractal Friar
  • *
  • Posts: 126
« Reply #10 on: January 03, 2019, 09:18:44 PM »
OK, from an old wip formula for UF (use Basic colouring from standard.ucl set to real with linear transfer density the same as the iteration count in the parameters i.e. detail - I have no recollection of why I called the formula "Two", I'm guessing it's based on similar info to that the author used in the images this thread refers to):

OK, I do not have UF. Can you make and post the image ?

Offline FractalDave

  • *
  • Uploader
  • *
  • Posts: 171
    • Makin Magic Fractals
« Reply #11 on: January 03, 2019, 11:19:19 PM »
OK, I do not have UF. Can you make and post the image ?

This is 2 layers, each just to the 21st iteration. One of every visible iteration boundary and one just rendering the boundary before "inside" and darkening the rest.
Only 21 iterations because here all point coordinate values are pre-calculated and stored before being added to the view.

It should be fairly easy to modify so that the colouring works fully on every boundary...

At some point between 2008 and now I did do a much more efficient inverse renderer that used some of my own speed-up ideas, but that was also before 2010 and it'll take me some time to find it !!!

Offline FractalDave

  • *
  • Uploader
  • *
  • Posts: 171
    • Makin Magic Fractals
« Reply #12 on: January 08, 2019, 01:48:48 AM »
Cracked it I think - pic attached.

Video (animating the start radius from 0.01 to 2.0) here: https://www.youtube.com/watch?v=yqeaWYf0zdE

Untidied code from test formula:

Code: [Select]
global:
  int i=0
  int j
  int k
  int l
  int m
  int n
  float o
  int s[@width+1,@height+1]
  complex v[2,round(2^@detail)];round(2^(24-@detail))]
  float vc[2,round(2^@detail)]
  int va[2,round(2^@detail)]
  int angles = round(2^@angles)
  float angles1 = angles
  int p[15]
  float q[282]
  float a
  float b
  float bb
  complex c
  complex zz
  b = 0
  repeat
;    if false;angles>0;282
      a = i*#pi*2.0/angles
;    else
;      a = q[i]
;    endif
    zz = v[0,0] = @start*(cos(a)+flip(sin(a)));cabs(sqrt(2-@const))*(cos(a)+flip(sin(a)));(0.5+sqrt(0.25+cabs(@const)))*(cos(a)+flip(sin(a)))
    vc[0,0]=vc[1,0]=vc[0,1]=vc[1,1]=0
    va[0,0]=round(i/angles)
;    if imag(c)<0
      vc[0,0]=vc[1,0]=i
;      print(i)
;    endif
    j = 0
    m = 0
    n = 1
    repeat
      k = 0
      l = round(2^j)
      ;a = 2^-(j+3)
      a = 1.0
      repeat
        c = v[n,l] = sqrt(v[m,k]-@const) ;1.0-(0.5-flip(0.5))*v[m,k,i];0.5*(sqrt(1+4.0*v[m,k,i])-1.0);
        o=vc[m,k]
        if o<0.5*angles
          o = o*0.5
        else
          o=0.5*(angles+o)
        endif
        if (o<0.5*angles && imag(c)<0) || (o>=0.5*angles && imag(c)>=0)
          o = angles-1-o
        endif
;        if o<0.5*angles
;          if imag(c)<0
;            o = 399/400;o + 0.5*angles
;          endif
;        elseif o>0.5*angles
;          if imag(c)>=0
;            o = o - 0.5*angles
;          endif
;        endif
        vc[n,l] = o
        if o*a>=1
          print(vc[m,k]," ",o," ",a)
        endif
        ;va[n,l]=b=0.5*va[n,l]+(1+atan2(c)/#pi)
        if @alliters || j==@detail-1
        if s[floor(0.5*@width+real(c)*@width*#magn/4.0),floor(0.5*@height-imag(c)*@height*#magn/3.0)-16]<=j
          s[floor(0.5*@width+real(c)*@width*#magn/4.0),floor(0.5*@height-imag(c)*@height*#magn/3.0)-16]=32*(angles*floor(o)+1)+j;(angles*2*round(b)+i+1)+j;(l-round(2^j))*angles+i+1
        endif
        endif
        c = v[n,k] = -c ;-1.0-(0.5+flip(0.5))*v[m,k,i];-0.5*(sqrt(1+4.0*v[m,k,i])+1.0);
        o=vc[m,k]
        if o<0.5*angles
          o = 0.5*(angles+o)
        else
          o = 0.5*o
        endif
        if o*a>=1
          print(vc[m,k]," ",o," ",a)
        endif
        if (o<0.5*angles && imag(c)<0) || (o>=0.5*angles && imag(c)>=0)
          o = angles-1-o
        endif
;        if o<0.5*angles
;          if imag(c)<0
;            o = o + 0.5*angles
;          endif
;        elseif o>0.5*angles
;          if imag(c)>=0
;            o = o - 0.5*angles
;          endif
;        endif
        vc[n,k] = o
;        if j<4
;        print(o)
;        endif
        ;va[n,k]=b=0.5*va[n,k]+1+atan2(c)/#pi
        if @alliters || j==@detail-1
        if s[floor(0.5*@width+real(c)*@width*#magn/4.0),floor(0.5*@height-imag(c)*@height*#magn/3.0)-16]<=j
          s[floor(0.5*@width+real(c)*@width*#magn/4.0),floor(0.5*@height-imag(c)*@height*#magn/3.0)-16]=32*(angles*floor(o)+1)+j;(angles*(2*round(b)+1)+i+1)+j;*round(2^(@detail-j-1))+i+1
        endif
        endif
        ;b = b + a
        l = l + 1
      until (k=k+1)>=round(2^j)
      k=0
      repeat
        print(i," ",vc[1,k]," ",vc[0,k])
      until (k=k+1)>=round(2^j)
      m = (m+1)%2
      n = (n+1)%2
    until (j=j+1)>=@detail
    ;print(i,"/",angles)
  until (i=i+1)>=angles
init:
  z = 0
  bool res = false
  float t
  if s[round(real(#screenpixel)),round(imag(#screenpixel))]>0
    res=true
    t = s[round(real(#screenpixel)),round(imag(#screenpixel))]
    t =(t-t%64)/32
    z = ((t%angles) + flip((t-t%angles)/angles))/angles;(angles*2^(@detail-1)))
  endif
loop:
bailout:
  !res
default:
  title = "Two"
  render=false
  complex param const
    default = (-0.12, 0.8)
    hint = "Julia seed"
  endparam
  bool param alliters
    default = false
  endparam
  float param start
    default = 2.0
  endparam
  int param detail
    default = 12
    max = 25
    min = 1
    hint = "max = 25 iterations"
  endparam
  int param angles
    caption = "#angles (2^)"
    default = 9
    min = 0
    max = 10
    hint = "max = 10 i.e. 2^10=1024"
  endparam
  int param width
    default = 640
    hint = "Should match image widrh"
  endparam
  int param height
    default = 480
    hint = "Should match image height"
  endparam
« Last Edit: January 08, 2019, 11:53:24 PM by FractalDave »

Offline claude

  • *
  • 3d
  • ****
  • Posts: 972
    • mathr.co.uk
« Reply #13 on: January 23, 2019, 09:25:53 AM »
not tested with disconnected Julia sets, it will probably explode (already crashes on 0+1i dentrite...)

algorithm is inverse iteration of a loop, do two passes per iteration (so you get double the number of points in the new loop, but then prune the list to avoid explosions).  picking the correct root to make the new loop smooth is almost trivially easy (assumes continuity).

Code: [Select]
// gcc -std=c99 -Wall -Wextra -pedantic -O3 -ffast-math -o julia-iim-rainbow julia-iim-rainbow.c -lm
// ./julia-iim-rainbow -0.125000000000000 0.649519052838329 > fat-rabbit.ppm

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define dpi 300
#define inches 6

#define size (dpi * inches)
unsigned char ppm[size][size][3];

static const double pi = 3.141592653589793;

double cnorm(double _Complex z)
{
  double x = creal(z);
  double y = cimag(z);
  return x * x + y * y;
}

inline void hsv2rgb(double hue, double sat, double bri, int *r, int *g, int *b)
{
  hue -= floor(hue);
  hue *= 6;
  int i = (int) floor(hue);
  double f = hue - i;
  if (! (i & 1))
    f = 1.0 - f;
  double m = bri * (1.0 - sat);
  double n = bri * (1.0 - sat * f);
  switch (i)
  {
  default:
  case 0: *r = 255 * bri;
    *g = 255 * n;
    *b = 255 * m;
    break;
  case 1: *r = 255 * n;
    *g = 255 * bri;
    *b = 255 * m;
    break;
  case 2: *r = 255 * m;
    *g = 255 * bri;
    *b = 255 * n;
    break;
  case 3: *r = 255 * m;
    *g = 255 * n;
    *b = 255 * bri;
    break;
  case 4: *r = 255 * n;
    *g = 255 * m;
    *b = 255 * bri;
    break;
  case 5: *r = 255 * bri;
    *g = 255 * m;
    *b = 255 * n;
  }
}

struct node
{
  struct node *next;
  double _Complex z;
  double h;
};

static struct node *freelist = 0;
struct node *allocate()
{
  if (freelist)
  {
    struct node *n = freelist;
    freelist = freelist->next;
    n->next = 0;
    n->z = 0;
    n->h = 0;
    return n;
  }
  else
  {
    return calloc(1, sizeof(struct node));
  }
}
void deallocate(struct node *n)
{
  n->next = freelist;
  freelist = n;
}

int main(int argc, char **argv)
{
  double _Complex c;
  if (argc > 2)
    c = atof(argv[1]) + I * atof(argv[2]);
  else if (argc > 1)
    c = atof(argv[1]);
  else
    c = 0;
  // repelling fixed point
  memset(&ppm[0][0][0], 0, size * size * 3);
  double scale = size / 5.0;
  double small = (0.25 * 0.25) / (scale * scale);
  double big = (1.0 * 1.0) / (scale * scale);
  // initialize loop
  struct node *loop = 0;
  for (int i = 0; i < 1 << 4; ++i)
  {
    double h = (i + 0.5) / (1 << 4);
    struct node *n = allocate();
    n->next = loop;
    n->z = 4 * cexp(2 * pi * I * h);
    n->h = h;
    loop = n;
  }
  // iterate
  for (int j = 0; j < 1 << 12; ++j)
  {
    // step loop twice (two branches)
    struct node *new_loop = 0;
    double _Complex z1 = 4;
    for (int pass = 0; pass < 2; ++pass)
    {
      struct node *old_loop = loop;
      while (old_loop)
      {
        // pick the root nearest the previous point in the new loop
        double _Complex z = csqrt(old_loop->z - c);
        double d1 = cnorm(z1 - z);
        double d2 = cnorm(z1 + z);
        struct node *n = allocate();
        n->next = new_loop;
        n->z = d1 < d2 ? z : -z;
        n->h = (old_loop->h + pass) * 0.5;
        new_loop = n;
        struct node *next = old_loop->next;
        if (pass == 1)
          deallocate(old_loop);
        old_loop = next;
        z1 = n->z;
      }
    }
    struct node *old_loop = new_loop;
    new_loop = 0;
    double _Complex z0 = old_loop->z;
    double _Complex h0 = old_loop->h;
    old_loop = old_loop->next;
    while (old_loop)
    {
      if (cnorm(old_loop->z - z0) < small)
      {
        // prune points that are too closely spaced
        struct node *next = old_loop->next;
        deallocate(old_loop);
        old_loop = next;
      }
      else if (cnorm(old_loop->z - z0) < big)
      {
        // move points that are reasonably spaced from old to new loop
        z0 = old_loop->z;
        h0 = old_loop->h;
        struct node *next = old_loop->next;
        old_loop->next = new_loop;
        new_loop = old_loop;
        old_loop = next;
      }
      else
      {
        // interpolate points that too distantly spaced
        struct node *next = allocate();
        next->next = old_loop;
        next->z = (z0 + old_loop->z) * 0.5;
        next->h = (h0 + old_loop->h) * 0.5;
        old_loop = next;
      }
    }
    loop = new_loop;
  }
  // plot loop
  while (loop)
  {
    int x = floor((creal(loop->z) + 2.5) * scale);
    int y = floor((cimag(loop->z) + 2.5) * scale);
    int r, g, b;
    hsv2rgb(loop->h, 1, 1, &r, &g, &b);
    ppm[y][x][0] = r;
    ppm[y][x][1] = g;
    ppm[y][x][2] = b;
    struct node *next = loop->next;
    deallocate(loop);
    loop = next;
  }
  fprintf(stdout, "P6\n# Julia set for %.18f + %.18f i\n%d %d\n255\n", creal(c), cimag(c), size, size);
  fwrite(&ppm[0][0][0], size * size * 3, 1, stdout);
  return 0;
}


xx
3D Julia

Started by trafassel on Fractal Image Gallery

0 Replies
68 Views
Last post August 22, 2018, 07:13:30 PM
by trafassel
question
A Julia set

Started by Matplotlib on Other

13 Replies
314 Views
Last post March 03, 2018, 06:09:26 AM
by SCORPION
xx
Julia In Motion

Started by utak3r on Fractal movie gallery

5 Replies
254 Views
Last post February 19, 2018, 01:56:11 AM
by utak3r
xx
Spiral Julia set in 3-D

Started by msltoe on Fractal Image Gallery

0 Replies
66 Views
Last post October 14, 2018, 07:08:15 AM
by msltoe
xx
Cantor Julia set

Started by Adam Majewski on Fractal Image Gallery

0 Replies
151 Views
Last post November 26, 2017, 07:19:35 PM
by Adam Majewski