Unveiling the Fractal Structure of Julia Sets with Lagrangian Descriptors

  • 60 Replies
  • 2727 Views

0 Members and 1 Guest are viewing this topic.

Online Adam Majewski

  • *
  • Fractal Fluff
  • *****
  • Posts: 396
« Reply #45 on: April 18, 2020, 11:14:00 AM »
Latest version. better but not perfect: boundaries of level sets in exterior are wide, in the interior not smooth. Help is welcome !!!

Code: [Select]
/*

  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero
 
 
 


 
  ==============================================
 
 
  Structure of a program or how to analyze the program
 
 
  ============== Image X ========================
 
  DrawImageOfX -> DrawPointOfX -> ComputeColorOfX
 
  first 2 functions are identical for every X
  check only last function =  ComputeColorOfX
  which computes color of one pixel !
 
 

   
  ==========================================

 
  ---------------------------------
  indent d.c
  default is gnu style
  -------------------



  c console progam
 
export  OMP_DISPLAY_ENV="TRUE"
  gcc d.c -lm -Wall -march=native -fopenmp
  time ./a.out > b.txt


  gcc d.c -lm -Wall -march=native -fopenmp


  time ./a.out

  time ./a.out >i.txt
  time ./a.out >e.txt
 
  =======================
  # gnuplot "i.plt"
set terminal svg enhanced background rgb 'white'
set xlabel "re(z)"
set ylabel "DLD"
set title "Relation between z and DLD in the interior of Julia set for c = -1"
set output "interior.svg"
plot "i.txt" with lines

  ----------------------
 
*/

#include <stdio.h>
#include <stdlib.h> // malloc
#include <string.h> // strcat
#include <math.h> // M_PI; needs -lm also
#include <complex.h>
#include <omp.h> // OpenMP

/* --------------------------------- global variables and consts ------------------------------------------------------------ */



// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1
//unsigned int ix, iy; // var
static unsigned int ixMin = 0; // Indexes of array starts from 0 not 1
static unsigned int ixMax; //
static unsigned int iWidth; // horizontal dimension of array

static unsigned int iyMin = 0; // Indexes of array starts from 0 not 1
static unsigned int iyMax; //

static unsigned int iHeight = 1000; // 
// The size of array has to be a positive constant integer
static unsigned int iSize; // = iWidth*iHeight;

// memmory 1D array
unsigned char *data;
//unsigned char *edge;
//unsigned char *edge2;

// unsigned int i; // var = index of 1D array
//static unsigned int iMin = 0; // Indexes of array starts from 0 not 1
static unsigned int iMax; // = i2Dsize-1  =
// The size of array has to be a positive constant integer
// unsigned int i1Dsize ; // = i2Dsize  = (iMax -iMin + 1) =  ;  1D array with the same size as 2D array


static const double ZxMin = -2.0; //-0.05;
static const double ZxMax =  2.0; //0.75;
static const double ZyMin = -2.0; //-0.1;
static const double ZyMax =  2.0; //0.7;
static double PixelWidth; // =(ZxMax-ZxMin)/ixMax;
static double PixelHeight; // =(ZyMax-ZyMin)/iyMax;
static double ratio;


// complex numbers of parametr plane
double complex c = -1.0; // parameter of function fc(z)=z^2 + c

double ER = 1e60;
double AR = 1e-60;



const int N = 1000; // fixed number : maximal number of iterations
double p  = 0.015625;  // 1/64






/* colors = shades of gray from 0 to 255 */
unsigned char iColorOfExterior = 250;
unsigned char iColorOfInterior = 200;
unsigned char iColorOfBoundary = 0;





/* ------------------------------------------ functions -------------------------------------------------------------*/





//------------------complex numbers -----------------------------------------------------





// from screen to world coordinate ; linear mapping
// uses global cons
double GiveZx ( int ix)
{
  return (ZxMin + ix * PixelWidth);
}

// uses globaal cons
double GiveZy (int iy) {
  return (ZyMax - iy * PixelHeight);
} // reverse y axis


complex double GiveZ( int ix, int iy){
  double Zx = GiveZx(ix);
  double Zy = GiveZy(iy);

  return Zx + Zy*I;




}




// ****************** DYNAMICS = trap tests ( target sets) ****************************


/* -----------  array functions = drawing -------------- */

/* gives position of 2D point (ix,iy) in 1D array  ; uses also global variable iWidth */
unsigned int Give_i (unsigned int ix, unsigned int iy)
{
  return ix + iy * iWidth;
}


// ***********************************************************************************************
// ********************** edge detection usung Sobel filter ***************************************
// ***************************************************************************************************

// from Source to Destination
int ComputeBoundaries(unsigned char S[], unsigned char D[])
{
 
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
  /* sobel filter */
  unsigned char G, Gh, Gv;
  // boundaries are in D  array ( global var )
 
  // clear D array
  memset(D, iColorOfExterior, iSize*sizeof(*D)); // for heap-allocated arrays, where N is the number of elements = FillArrayWithColor(D , iColorOfExterior);
 
  // printf(" find boundaries in S array using  Sobel filter\n");   
#pragma omp parallel for schedule(dynamic) private(i,iY,iX,Gv,Gh,G) shared(iyMax,ixMax)
  for(iY=1;iY<iyMax-1;++iY){
    for(iX=1;iX<ixMax-1;++iX){
      Gv= S[Give_i(iX-1,iY+1)] + 2*S[Give_i(iX,iY+1)] + S[Give_i(iX-1,iY+1)] - S[Give_i(iX-1,iY-1)] - 2*S[Give_i(iX-1,iY)] - S[Give_i(iX+1,iY-1)];
      Gh= S[Give_i(iX+1,iY+1)] + 2*S[Give_i(iX+1,iY)] + S[Give_i(iX-1,iY-1)] - S[Give_i(iX+1,iY-1)] - 2*S[Give_i(iX-1,iY)] - S[Give_i(iX-1,iY-1)];
      G = sqrt(Gh*Gh + Gv*Gv);
      i= Give_i(iX,iY); /* compute index of 1D array from indices of 2D array */
      if (G==0) {D[i]=255;} /* background */
      else {D[i]=0;}  /* boundary */
    }
  }
 
   
 
  return 0;
}



// copy from Source to Destination
int CopyBoundaries(unsigned char S[],  unsigned char D[])
{
 
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
 
 
  //printf("copy boundaries from S array to D array \n");
  for(iY=1;iY<iyMax-1;++iY)
    for(iX=1;iX<ixMax-1;++iX)
      {i= Give_i(iX,iY); if (S[i]==0) D[i]=0;}
 
 
 
  return 0;
}







// ***************************************************************************************************************************
// ************************** DLD/J*****************************************
// ****************************************************************************************************************************



/* partial pnorm
   input: z , zn = f(z), p
   output ppn
   
   
*/
double ppnorm( complex double z, complex double zn, double p){

double s[2][3]; // array for 2 points on the Riemann sphere
int j;
double d; // denominator
double x;
double y;

double ds;
double ppn = 0.0;

// map from complex plane to riemann sphere
// z
x = creal(z);
y = cimag(z);
d = x*x + y*y + 1.0;

s[0][0] = (2.0*x)/d;
s[0][1] = (2.0*y)/d; 
s[0][2] = (d-2.0)/d; // (x^2 + y^2 - 1)/d

// zn
x = creal(zn);
y = cimag(zn);
d = x*x + y*y + 1.0;
s[1][0] = (2.0*x)/d;
s[1][1] = (2.0*y)/d; 
s[1][2] = (d-2.0)/d; // (x^2 + y^2 - 1)/d

// sum
for (j=0; j <3; ++j){
ds = fabs(s[1][j] - s[0][j]);
//  normal:  neither zero, subnormal, infinite, nor NaN
//if (fpclassify (ds) !=FP_INFINITE)
//if (isnormal(ds))
// it is solved by if (cabs(z) > 1e60 ) break; procedure in parent function
ppn += pow(ds,p); // |ds|^p
// else {ppn = 10000.0; printf("ds = infty\t");} //

}


return ppn;







}

// DLD = Discret Lagrangian Descriptior
double lagrangian( complex double z0, complex double c, int iMax, double p ){

int i; // number of iteration
double d = 0.0; // DLD = sum
double ppn; // partial pnorm
complex double z = z0;
complex double zn; // next z


if (cabs(z) < AR || cabs(z +1)< AR) return 5.0; // for z= 0.0 d = inf


for (i=0; i<iMax; ++i){




zn = z*z +c; // complex iteration
ppn = ppnorm(z, zn, p);
d += ppn; // sum
//
z = zn;

//if (! isnormal(d)) { return 0.0; } // not works
if (cabs(z) > ER ) break; // exterior : big values produces NAN error in ppnorm computing
if (cabs(z) < AR || cabs(z +1)< AR)
{ // interior
d = -d;
break;

}


}


d =  d/((double)i); // averaging not summation
if (d<0.0) {// interior
d = 2.5 - d;


}
//d = d - (int)d; // fractional part
return d;




}





unsigned char ComputeColorOfDLD(complex double z){

 
  //double cabsz;
  int iColor;
  double d;
 

  d = lagrangian(z,c, N,p);
 
 
    iColor = (int)(d*255)  % 255; // nMax or lower walues in denominator
 
 
  return (unsigned char) iColor;


}



// plots raster point (ix,iy)
int DrawPointOfDLD (unsigned char A[], int ix, int iy)
{
  int i; /* index of 1D array */
  unsigned char iColor;
  complex double z;


  i = Give_i (ix, iy); /* compute index of 1D array from indices of 2D array */
  z = GiveZ(ix,iy);
  iColor = ComputeColorOfDLD(z);
  A[i] = iColor ; // interior
 
  return 0;
}




// fill array
// uses global var :  ...
// scanning complex plane
int DrawImagerOfDLD (unsigned char A[])
{
  unsigned int ix, iy; // pixel coordinate

  //printf("compute image \n");
  // for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  for (iy = iyMin; iy <= iyMax; ++iy){
    //printf (" %d from %d \r", iy, iyMax); //info
    for (ix = ixMin; ix <= ixMax; ++ix)
      DrawPointOfDLD(A, ix, iy); // 
  }

  return 0;
}





// test how values changes to tune color
int test_interior(){

complex double z = 0.0;
complex double dz = 0.001;



printf("# z d\n"); // gnuplot
while (cabs(z) < 0.6){

double d = lagrangian(z, c, N, p);
int iColor = ComputeColorOfDLD(z);

// printf(" z = %.16f d = %.16f color = %d \n",creal(z), d, iColor);
printf(" %.16f %.16f \n",creal(z), d); // gnuplot
z += dz;
}

//
double d0 = lagrangian(0.00000, c, N, p);
double db = lagrangian(0.63125, c, N, p);
double dd = d0 - db;
printf("d0 - db  = %.16f - %.16f = %.16f\n",d0, db, dd);

return 0;


}
 
 
 

// test how values changes to tune color
int test_exterior(){

complex double z;
complex double z0 = 1.62;
complex double z1 = 3.0;
complex double dz = 0.001;


z = z0;
printf("# z d\n"); // gnuplot
while (creal(z) < creal(z1)){

double d = lagrangian(z, c, N, p);
int iColor = ComputeColorOfDLD(z);

// printf(" z = %.16f d = %.16f color = %d \n",creal(z), d, iColor);
printf(" %.16f %.16f \n",creal(z), d); // gnuplot
z += dz;
}

//
double d0 = lagrangian(z0, c, N, p);
double d1 = lagrangian(z1, c, N, p);
double dd = d0 - d1;
printf("d0 - d1  = %.16f - %.16f = %.16f\n",d0, d1, dd);

return 0;


}
 
 








// *******************************************************************************************
// ********************************** save A array to pgm file ****************************
// *********************************************************************************************

int SaveArray2PGMFile( unsigned char A[], double k, char* comment )
{
 
  FILE * fp;
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name [100]; /* name of file */
  snprintf(name, sizeof name, "%.0f", k); /*  */
  char *filename =strncat(name,".pgm", 4);
 
 
 
  // save image to the pgm file
  fp= fopen(filename,"wb"); // create new file,give it a name and open it in binary mode
  fprintf(fp,"P5\n # %s\n %u %u\n %u\n", comment, iWidth, iHeight, MaxColorComponentValue);  // write header to the file
  fwrite(A,iSize,1,fp);  // write array with image data bytes to the file in one step
  fclose(fp);
 
  // info
  printf("File %s saved ", filename);
  if (comment == NULL || strlen(comment) ==0) 
    printf("\n");
  else printf (". Comment = %s \n", comment);

  return 0;
}







int PrintInfoAboutProgam()
{

 
  // display info messages
  printf ("Numerical approximation of Julia set for fc(z)= z^2 + c \n");
  //printf ("iPeriodParent = %d \n", iPeriodParent);
  //printf ("iPeriodOfChild  = %d \n", iPeriodChild);
  printf ("parameter c = ( %.16f ; %.16f ) \n", creal(c), cimag(c));
 
  printf ("Image Width = %f in world coordinate\n", ZxMax - ZxMin);
  printf ("PixelWidth = %f \n", PixelWidth);
 
 
  // image corners in world coordinate
  // center and radius
  // center and zoom
  // GradientRepetition
  printf ("Maximal number of iterations = iterMax = %d \n", N);
  printf ("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  //
  printf("gcc version: %d.%d.%d\n",__GNUC__,__GNUC_MINOR__,__GNUC_PATCHLEVEL__); // https://stackoverflow.com/questions/20389193/how-do-i-check-my-gcc-c-compiler-version-for-my-eclipse
  // OpenMP version is diplayed in the console
  return 0;
}






// *****************************************************************************
//;;;;;;;;;;;;;;;;;;;;;;  setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
// **************************************************************************************

int setup ()
{

  printf ("setup start\n");
   
 
 
 
 

  /* 2D array ranges */
 
  iWidth = iHeight;
  iSize = iWidth * iHeight; // size = number of points in array
  // iy
  iyMax = iHeight - 1; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].
  //ix

  ixMax = iWidth - 1;

  /* 1D array ranges */
  // i1Dsize = i2Dsize; // 1D array with the same size as 2D array
  iMax = iSize - 1; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].

  /* Pixel sizes */
  PixelWidth = (ZxMax - ZxMin) / ixMax; //  ixMax = (iWidth-1)  step between pixels in world coordinate
  PixelHeight = (ZyMax - ZyMin) / iyMax;
  ratio = ((ZxMax - ZxMin) / (ZyMax - ZyMin)) / ((double) iWidth / (double) iHeight); // it should be 1.000 ...

   

 
   
  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc (iSize * sizeof (unsigned char));
   
 
  if (data == NULL  ){
    fprintf (stderr, " Could not allocate memory");
    return 1;
  }

 
 
 
 
 
 
  printf (" end of setup \n");

  return 0;

} // ;;;;;;;;;;;;;;;;;;;;;;;;; end of the setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;




int end(){


  printf (" allways free memory (deallocate )  to avoid memory leaks \n"); // https://en.wikipedia.org/wiki/C_dynamic_memory_allocation
  free (data);
 
 
  PrintInfoAboutProgam();
  return 0;

}

// ********************************************************************************************************************
/* -----------------------------------------  main   -------------------------------------------------------------*/
// ********************************************************************************************************************

int main () {
  setup ();
 
 
 
  DrawImagerOfDLD(data);
  SaveArray2PGMFile (data, iWidth, "DLD/J");


  //test_exterior();
  test_interior();
 
  end();

  return 0;
}

« Last Edit: April 18, 2020, 07:03:20 PM by Adam Majewski »

Online Adam Majewski

  • *
  • Fractal Fluff
  • *****
  • Posts: 396
« Reply #46 on: April 19, 2020, 10:08:05 AM »
... averaging keeps the coloring stable if maxiters is changed but can lead to low variation over the image



Is it related ?

Offline pauldelbrot

  • *
  • 3f
  • ******
  • Posts: 2295
« Reply #47 on: April 21, 2020, 01:11:13 AM »
More supernova parameter space; pure DLD coloring. A(0) and A(∞) behave so alike, the DLD coloring treats them identically, but they alternate, switching when one crosses one of the maroon fibers. The finite attractors, when present, behave differently enough that the DLD coloring assigns from a different region of the gradient, producing blue to turquoise minibrots and associated features.

Online Adam Majewski

  • *
  • Fractal Fluff
  • *****
  • Posts: 396
« Reply #48 on: May 01, 2020, 08:32:44 PM »

Online Adam Majewski

  • *
  • Fractal Fluff
  • *****
  • Posts: 396
« Reply #49 on: May 01, 2020, 10:12:12 PM »
« Last Edit: May 03, 2020, 04:24:32 PM by Adam Majewski »

Offline C0ryMcG

  • *
  • Uploader
  • *
  • Posts: 241
« Reply #50 on: May 03, 2020, 08:35:47 AM »
I think I've finally got this working, or at least I'm very close to it...

Online Adam Majewski

  • *
  • Fractal Fluff
  • *****
  • Posts: 396
« Reply #51 on: May 03, 2020, 04:26:09 PM »
cool.
I like the exterior where field lines and equipotential curves ( boundaries of level sets) are seen

Can you made similar effect for interior?

Offline C0ryMcG

  • *
  • Uploader
  • *
  • Posts: 241
« Reply #52 on: May 04, 2020, 06:19:20 AM »
 
cool.
I like the exterior where field lines and equipotential curves ( boundaries of level sets) are seen

Can you made similar effect for interior?

I haven't made that to work yet, although an idea I have is to inverse the orbit locations when the orbit is shown to not escape and run the same calculations on that inversion.
Meanwhile, here is an example of an inverted fractal with this coloring, so the details end up in the center.. I like the result.

Online Adam Majewski

  • *
  • Fractal Fluff
  • *****
  • Posts: 396
« Reply #53 on: May 04, 2020, 09:49:46 PM »
Fat basilica - field lines (interior)

Here I see some error: all curves shoud start and end at parabolic fixed point ant it's preimages .


Offline C0ryMcG

  • *
  • Uploader
  • *
  • Posts: 241
« Reply #54 on: May 16, 2020, 04:25:10 AM »
An update... For reasons that are not entirely clear to me, I noticed that if I changed the Power of my mandelbrot equation to anything that was not an integer multiple of 2, the equipotential lines in my renders disappeared, and the color contrast dropped. I could use z^2.001 and z^1.999 and get clean (but boring) renders looking almost like regular escape time coloring. But z^2 or z^4 would get those fantastic stripes seen in my previous images.

This has something to do with my code's error avoidance. It won't add to the orbit array if the result of the iteration is NaN.  And I was using the array length as the average formula's divisor... I changed this divisor to the Iteration Count and my equipotential lines disappeared for all numbers, including 2 and multiples. On the one hand, this means my calculated colors are more accurate for the formula.
On the other, my images are less exciting.  So I'm a bit torn. I'd appreciate if anyone has any insight to why Javascript might not like certain escaping points when the power is 2, but has no trouble with any nearby power. Maybe if I could understand it, I could control it and have customizable lines.
Does anyone have an idea of why Javascript mi

Offline marcm200

  • *
  • 3c
  • ***
  • Posts: 855
« Reply #55 on: May 17, 2020, 09:23:47 AM »
Sorry, I don't have an answer, but some questions (it's like another detective story :) ): How did you implement the iteration formula for the integer exponents? As multiplications - and for rational powers as exp/log? Or is everything exp/log? In case, the former is true, what does your program output if you use exp/log and then 2.0 as an exponent? What happens with odd integers?

Offline C0ryMcG

  • *
  • Uploader
  • *
  • Posts: 241
« Reply #56 on: May 18, 2020, 12:12:04 AM »
Ah, you provide some good hints... As I think you suspected, my Mandelbrot iteration formula treats z^2 (and z^4) differently than z^p when p!=2 or 4. This was done to speed things up for the default rendering, since I have a much simpler formula for complexSquare than complexPower, and ^4 can be done as Square(Square()). I did not test power of 6 until now, and assumed it would behave like 2 or 4, which was incorrect.

So, then, the difference between the two functions... Mathmatically they are the same. But after some testing, I saw that the one function returns NaN when Z gets out of bounds of the float, while the other returns Infinity or -Infinity. One of these gets added to the array, the other does not. So if the power is 2 or 4, the array of resulting points in the orbit returns shortened. If it's anything else, it returns at full length, with a lot of 'infinities' added on to the end.

So thanks for the nudge in the right direction! Now I'm going to go back to facing the aesthetic side of the problem, which is often so much harder to solve..

Offline QBongo

  • *
  • Fractal Freshman
  • *
  • Posts: 1
« Reply #57 on: June 13, 2020, 01:03:41 AM »
hi, can you tell me what language this is?

Code: [Select]
DLD {
; Based on https://arxiv.org/pdf/2001.08937.pdf
init:
  float sum = 0.0
  float lastx = 0.0
  float lasty = 0.0
  float lastz = 0.0
  int i = 0
loop:
  float d = |#z|
  float dd = 1/(d + 1)
  float xx = 2*real(#z)*dd
  float yy = 2*imag(#z)*dd
  float zz = (d - 1)*dd     ; Riemann sphere coordinates
  IF (i > 0)
    sum = sum + (xx - lastx)^@power + (yy - lasty)^@power + (zz - lastz)^@power
  ENDIF
  i = i + 1
  lastx = xx
  lasty = yy
  lastz = zz
final:
  #index = sum/(i - 1)
default:
  title = "Discrete Langrangian Descriptors"
  param power
    caption = "Power"
    default = 0.25
    hint = "Affects the steepness of the gradient near singularities (fractal features like a Julia set)"
    min = 0.0
  endparam
}

Offline pauldelbrot

  • *
  • 3f
  • ******
  • Posts: 2295
« Reply #58 on: June 13, 2020, 02:09:16 AM »
Ultrafractal's formula language. Will work in UF6 and probably works in UF5, 4.

Offline 3DickUlus

  • *
  • Administrator
  • *******
  • Posts: 1866
    • Digilantism
« Reply #59 on: June 22, 2020, 10:58:00 PM »
Thanks for sharing Adam  :thumbs:

field lines and equipotential curves ? maybe so tight they disappear?

Fragment code at https://fractalforums.org/code-snippets-fragments/74/lagrangian-descriptors-fragment-code/3612/msg22392
« Last Edit: June 23, 2020, 04:35:47 AM by 3DickUlus »


clip
Lagrangian Descriptors Fragment code

Started by 3DickUlus on Code Snippets (fragments)

10 Replies
439 Views
Last post June 30, 2020, 01:32:59 AM
by 3DickUlus
xx
Embedded Julia sets

Started by pauldelbrot on Image Threads

108 Replies
6878 Views
Last post November 16, 2020, 02:26:14 AM
by pauldelbrot
clip
Parabolic Julia sets

Started by marcm200 on Fractal Mathematics And New Theories

82 Replies
3400 Views
Last post November 11, 2020, 02:29:07 PM
by marcm200
clip
Quadratic Julia sets

Started by pauldelbrot on Image Threads

26 Replies
1381 Views
Last post Today at 01:40:17 PM
by pauldelbrot
clip
3d Julia sets: True shape

Started by marcm200 on Image Threads

31 Replies
1325 Views
Last post January 20, 2020, 11:42:23 AM
by marcm200