Logo by KRAFTWERK - Contribute your own Logo!

END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

it was a great time but no longer maintainable by c.Kleinhuis contact him for any data retrieval,
thanks and see you perhaps in 10 years again

this forum will stay online for reference
News: Check out the originating "3d Mandelbulb" thread here
 
*
Welcome, Guest. Please login or register. March 29, 2024, 11:35:21 AM


Login with username, password and session length


The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!


Pages: [1]   Go Down
  Print  
Share this topic on DiggShare this topic on FacebookShare this topic on GoogleShare this topic on RedditShare this topic on StumbleUponShare this topic on Twitter
Author Topic: DEM/M in c and shadertoy  (Read 6398 times)
0 Members and 1 Guest are viewing this topic.
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« on: January 11, 2014, 01:21:15 AM »

Hi,

There is a nice Mandelbrot set in Shadertoy : Mandelbrot - distance
http://www.shadertoy.com/view/lsX3W4

with code by inigo quilez 

Code:
void main(void)
{
    vec2 p = -1.0 + 2.0 * gl_FragCoord.xy / iResolution.xy;
    p.x *= iResolution.x/iResolution.y;

    // animation
float tz = 0.5 - 0.5*cos(0.225*iGlobalTime);
    float zoo = pow( 0.5, 13.0*tz );
vec2 cc = vec2(-0.05,.6805) + p*zoo;

    // iterate
    vec2 z  = vec2(0.0);
    float m2 = 0.0;
    float co = 0.0;
    vec2 dz = vec2(0.0);
    for( int i=0; i<256; i++ )
    {
        if( m2>1024.0 ) continue;

// Z' -> 2·Z·Z' + 1
        dz = 2.0*vec2(z.x*dz.x-z.y*dz.y, z.x*dz.y + z.y*dz.x ) + vec2(1.0,0.0);

        // Z -> Z² + c
        z = cc + vec2( z.x*z.x - z.y*z.y, 2.0*z.x*z.y );

        m2 = dot(z,z);
        co += 1.0;
    }

    // distance
// d(c) = |Z|·log|Z|/|Z'|
float d = 0.0;
    if( co<256.0 ) d = sqrt( dot(z,z)/dot(dz,dz) )*log(dot(z,z));


    // do some soft coloring based on distance
d = clamp( 4.0*d/zoo, 0.0, 1.0 );
d = pow( d, 0.25 );
    vec3 col = vec3( d );
   
    gl_FragColor = vec4( col, 1.0 );
}

I try to do similar ( without animation ) in c  ( code below ) but the effect is not good .
How can I change it ?

TIA

 c code below :
Code:
// gcc m.c  -lm -Wall -march=native
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>



/* iXmax/iYmax = 1 */
unsigned int iSide = 1000; /* side of rectangle in pixels */
unsigned int iXmax; // ((int)(m*iSide)) /* height of image in pixels */
unsigned int iYmax ; //= iSide;
unsigned int iLength; // = (iXmax*iYmax) /* number of pixels */
/*  world ( double) coordinate */
double dSide =  1.5;
double CxMin ; // = 0.0;
double CxMax ; // =(m*dSide);
double CyMin ; //= dSide;
double CyMax ; // = dSide;
/* (CxMax-CxMin)/(CyMax-CyMin)==iXmax/iYmax = 1 */

unsigned int IterationMax; // = (iXmax*100) /* proportional to resolution of picture */

double PixelWidth; //=  ((CxMax-CxMin)/iXmax)
double PixelHeight ;//= ((CyMax-CyMin)/iYmax)

double CDistanceMax ; //= PixelWidth; /* proportional to pixel size */


/* fc(z) = z*z + c */



double EscapeRadius = 33.0;  /* radius of circle around origin; its complement is a target set for escaping points */
double ER2 ; //= (EscapeRadius*EscapeRadius)


/* colors = shades of gray from 0=black  to 255=white */
unsigned char iExterior = 245; /* exterior of Julia set */
unsigned char iBoundary = 0; /* border , boundary*/
unsigned char iInterior = 180;







unsigned int f(unsigned int _iX, unsigned int _iY)
/*
   gives position of point (iX,iY) in 1D array  ; uses also global variables
   it does not check if index is good  so memory error is possible
*/
{return (_iX + (iYmax-_iY-1)*iXmax );}








/*
 estimates distance from point c to nearest point in Julia  set
 for Fc(z)= z*z + c
 z(n+1) = Fc(zn) 
 this function is based on function  mndlbrot::dist  from  mndlbrot.cpp
 from program mandel by Wolf Jung (GNU GPL )
 http://www.mndynamics.com/indexp.html

Hyunsuk Kim  :
For Julia sets, z is the variable and c is a constant. Therefore df[n+1](z)/dz = 2*f[n]*f'[n] -- you don't add 1.

For the Mandelbrot set on the parameter plane, you start at z=0 and c becomes the variable. df[n+1](c)/dc = 2*f[n]*f'[n] + 1.
http://iquilezles.org/www/articles/distancefractals/distancefractals.htm

 */

// boolean Escape time and DEM/M in one loop 
unsigned char GiveColor( double Cx, double Cy, int iMax, double DistanceMax)
{

  // C = Cx + Cy* I = point of parameter c-plane

  int i; /* iteration number */

  double Zx, Zy; /* Z = Zx + Zy*I  point of dynamical plane */
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  double temp;
  double absZ2;

  // http://en.wikipedia.org/wiki/Complex_quadratic_polynomial
  // first derivative  of fc(zcr) with respect to c =  dZ = dfc(zcr)/dc = 2*Z*dZ  = dZx + dZy*I
  double dZx = 0.0;
  double dZy = 0.0;
  double absdZ2; // = abs(dZ)* abs(dZ) = dZx*dZx + dZy*dZy

  double distance;
  unsigned char color = iInterior;

 

  /* initial value of orbit  = critical point zcr = 0 */
  Zx=0.0;
  Zy=0.0;
  //
  Zx2 = Zx*Zx;
  Zy2 = Zy*Zy;
  absZ2= Zx*Zx + Zy*Zy;
  absdZ2= dZx*dZx + dZy*dZy;
  // iteration of critical point z= 0 on the dynamical z-plane       
  for (i=0; i<iMax; i++)
    { // check if not escaping : abs(z)>ER
      if (absZ2  > ER2 ) { color=  iExterior; break ; } // exterior when escapes
      //if (absdZ2 > 1e60) { return iInterior; } //  interior when derivative explodes
      // z = fc(z) = z*z + c
      Zy=2*Zx*Zy + Cy;
      Zx=Zx2-Zy2 +Cx;
     
     /* first derivative   zp = 2*z*zp  = xp + yp*i; */
      temp = 2*(Zx*dZx - Zy*dZy) + 1.0 ; /*  */
      dZy = 2*(Zx*dZy + Zy*dZx);
      dZx = temp;

      // abs
      Zx2 = Zx*Zx;     
      Zy2 = Zy*Zy;
      absZ2= Zx2 + Zy2; // nz = x*x + y*y;
      absdZ2= dZx*dZx + dZy*dZy; //
    };

  if (i<iMax)
   {
     
     distance = sqrt(absZ2/absdZ2)*log(absZ2); //
     distance = pow( distance, 0.25 );
     //if (distance > 0.0 && distance<DistanceMax*)
       //   { //color = ((int)(255.0*distance)) % 255;

            // d = sqrt( dot(z,z)/dot(dz,dz) )*log(dot(z,z));
        // do some soft coloring based on distance
    //d = pow( 4.0*d/zoom, 0.25 );
            //if (distance < (DistanceMax)) color = iBoundary ; // bondary = black
            // near boundary = shades of gray
            color = 255-255*(int)(distance/DistanceMax);
        //}   
    }
   // if (nz < bailout) return 1; //still not escaping after iteration , rare
  // if (absdZ2 < absZ2) color= iExterior;  //includes escaping through 0 // som eiterior points are coded as exterior = error
   
  return color;
} //0.000007
 


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

int main(){

  unsigned int iX,iY, /* indices of 2D virtual array (image) = integer coordinate */
    i; /* index of 1D array  */
  double Cx, Cy;
 
  //int  ExtLastIteration; //
  unsigned char color;


  printf(" Setup \n");
  iXmax= iSide; /* height of image in pixels */
  iYmax  = iSide;
  iLength = iXmax*iYmax; /* number of pixels */
  CxMin = -2.25;
  CxMax = 0.75  ;
  CyMin = -dSide;
  CyMax = dSide;
  IterationMax = 1000 ;
 PixelWidth=  (CxMax-CxMin)/iXmax;
 PixelHeight = (CyMax-CyMin)/iYmax;
  ER2 = EscapeRadius*EscapeRadius;

  /* dynamic 1D array for colors ( shades of gray ) */
  unsigned char *data;
  data = malloc( iLength * sizeof(unsigned char) );
  if (data == NULL )
    {
      fprintf(stderr," Could not allocate memory");
      return 1;
    }
   
  printf(" compute color \n");
  for(iY=0;iY<iYmax;++iY){
    Cy=CyMin + iY*PixelHeight; /*  */
    printf("row %u from %u \n",iY, iYmax);   
    for(iX=0;iX<iXmax;++iX){
      Cx=CxMin + iX*PixelWidth;
      color=GiveColor(Cx, Cy, IterationMax, PixelWidth);
     
      i= f(iX,iY); /* compute index of 1D array from indices of 2D array */
      data[i]=color;  /* change the color */
     
    }
  }


 

   

  printf(" save  data array to the pgm file \n");
  unsigned int length = 30;
  char filename[length] ;
  snprintf(filename, length, "%.0fe%u%s", EscapeRadius,IterationMax, ".pgm");
  char *comment="#  ";/* comment should start with # */
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  /* save image to pgm file  */
  FILE * fp;     
  fp = fopen(filename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"P5\n %s\n %u\n %u\n %u\n",comment,iXmax,iYmax,MaxColorComponentValue);  /*write header to the file*/
  fwrite(data,iLength,1,fp);  /*write image data bytes to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);

  /* ---------- text file  -------------------------------------*/
  char tfilename[length] ;
  printf(" save  graphic data to the text file \n");
  snprintf(tfilename, length, "%.0fe%u%s",EscapeRadius, iXmax, ".txt");
   /* save info  to the file  */     
  fp = fopen(tfilename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"IterationMax = %d \n", IterationMax);
  fprintf(fp,"EscapeRadius = %f ER2 = %f \n", EscapeRadius, ER2);
 

  fprintf(fp,"\n" );  /* */
  fprintf(fp,"C plane : \n" );  /* */
  fprintf(fp,"dSide = %f \n", dSide);
  fprintf(fp,"PixelWidth = %f ; PixelHeight = %f \n", PixelHeight, PixelWidth);
  fprintf(fp," CxMin = %f \n", CxMin);  /* */
  fprintf(fp," CxMax = %f \n", CxMax);  /* */
  fprintf(fp," CyMin = %f \n", CyMin);  /* */
  fprintf(fp," CyMax = %f \n", CyMax);  /* */
  fprintf(fp,"\n" );
  //
  printf("File %s saved. \n", tfilename);
  fclose(fp);


  printf(" allways free memory \n ");
  free(data);
 
 

  return 0;
}
Logged
xenodreambuie
Conqueror
*******
Posts: 124



WWW
« Reply #1 on: January 11, 2014, 07:46:23 AM »

I try to do similar ( without animation ) in c  ( code below ) but the effect is not good .
How can I change it ?

One obvious problem is that in the loop, the derivative should be calculated before the new z.
Logged

Regards, Garth
http://xenodream.com
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« Reply #2 on: January 11, 2014, 09:02:02 AM »


  • One obvious problem is that in the loop, the derivative should be calculated before the new z.
    OK. It looks better now.
  • In GSLS code color is 32 bit and every component is in range from 0.0 to 1.0. In c code color is 8 bit = 1 component and is in range  0 to 255. SO I have to mulitply  : color = 255-255*(int)(distance);

Still not the same

Help is wellcome

Code:
/*

 

  to compile :
  gcc m.c  -lm -Wall -march=native
  to run ( Linux console) :
  time ./a.out


*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>



/* iXmax/iYmax = 1 */
unsigned int iSide = 1000; /* side of rectangle in pixels */
unsigned int iXmax; // ((int)(m*iSide)) /* height of image in pixels */
unsigned int iYmax ; //= iSide;
unsigned int iLength; // = (iXmax*iYmax) /* number of pixels */
/*  world ( double) coordinate */
double dSide =  1.5;
double CxMin ; // = 0.0;
double CxMax ; // =(m*dSide);
double CyMin ; //= dSide;
double CyMax ; // = dSide;
/* (CxMax-CxMin)/(CyMax-CyMin)==iXmax/iYmax = 1 */

unsigned int IterationMax; // = (iXmax*100) /* proportional to resolution of picture */

double PixelWidth; //=  ((CxMax-CxMin)/iXmax)
double PixelHeight ;//= ((CyMax-CyMin)/iYmax)

double CDistanceMax ; //= PixelWidth; /* proportional to pixel size */


/* fc(z) = z*z + c */



double EscapeRadius = 33.0;  /* radius of circle around origin; its complement is a target set for escaping points */
double ER2 ; //= (EscapeRadius*EscapeRadius)


/* colors = shades of gray from 0=black  to 255=white */
unsigned char iExterior = 245; /* exterior of Julia set */
unsigned char iBoundary = 0; /* border , boundary*/
unsigned char iInterior = 180;







unsigned int f(unsigned int _iX, unsigned int _iY)
/*
   gives position of point (iX,iY) in 1D array  ; uses also global variables
   it does not check if index is good  so memory error is possible
*/
{return (_iX + (iYmax-_iY-1)*iXmax );}








/*
 estimates distance from point c to nearest point in Julia  set
 for Fc(z)= z*z + c
 z(n+1) = Fc(zn) 
 this function is based on function  mndlbrot::dist  from  mndlbrot.cpp
 from program mandel by Wolf Jung (GNU GPL )
 http://www.mndynamics.com/indexp.html

Hyunsuk Kim  :
For Julia sets, z is the variable and c is a constant. Therefore df[n+1](z)/dz = 2*f[n]*f'[n] -- you don't add 1.

For the Mandelbrot set on the parameter plane, you start at z=0 and c becomes the variable. df[n+1](c)/dc = 2*f[n]*f'[n] + 1.
http://iquilezles.org/www/articles/distancefractals/distancefractals.htm

 */

// boolean Escape time and DEM/M in one loop 
unsigned char GiveColor( double Cx, double Cy, int iMax, double DistanceMax)
{

  // C = Cx + Cy* I = point of parameter c-plane

  int i; /* iteration number */

  double Zx, Zy; /* Z = Zx + Zy*I  point of dynamical plane */
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  double temp;
  double absZ2;

  // http://en.wikipedia.org/wiki/Complex_quadratic_polynomial
  // first derivative  of fc(zcr) with respect to c =  dZ = dfc(zcr)/dc = 2*Z*dZ  = dZx + dZy*I
  double dZx = 0.0;
  double dZy = 0.0;
  double absdZ2; // = abs(dZ)* abs(dZ) = dZx*dZx + dZy*dZy

  double distance;
  unsigned char color = iInterior;

 

  /* initial value of orbit  = critical point zcr = 0 */
  Zx=0.0;
  Zy=0.0;
  //
  Zx2 = Zx*Zx;
  Zy2 = Zy*Zy;
  absZ2= Zx*Zx + Zy*Zy;
  absdZ2= dZx*dZx + dZy*dZy;
  // iteration of critical point z= 0 on the dynamical z-plane       
  for (i=0; i<iMax; i++)
    { // check if not escaping : abs(z)>ER
      if (absZ2  > ER2 ) { color=  iExterior; break ; } // exterior when escapes
      //if (absdZ2 > 1e60) { return iInterior; } //  interior when derivative explodes
     

       // in the loop, the derivative should be calculated before the new z
      /* first derivative   zp = 2*z*zp  = xp + yp*i; */
      temp = 2*(Zx*dZx - Zy*dZy) + 1.0 ; /*  */
      dZy = 2*(Zx*dZy + Zy*dZx);
      dZx = temp;

      // z = fc(z) = z*z + c
      Zy=2*Zx*Zy + Cy;
      Zx=Zx2-Zy2 +Cx;
     
     

      // abs
      Zx2 = Zx*Zx;     
      Zy2 = Zy*Zy;
      absZ2= Zx2 + Zy2; // nz = x*x + y*y;
      absdZ2= dZx*dZx + dZy*dZy; //
    };

  if (i<iMax)
   {
     
     distance = sqrt(absZ2/absdZ2)*log(absZ2); //
     distance = pow( 4.0*distance, 0.25 );
     //if (distance > 0.0 && distance<DistanceMax*)
       //   { //color = ((int)(255.0*distance)) % 255;

            // d = sqrt( dot(z,z)/dot(dz,dz) )*log(dot(z,z));
        // do some soft coloring based on distance
    //d = pow( 4.0*d/zoom, 0.25 );
            //if (distance < (DistanceMax)) color = iBoundary ; // bondary = black
            // near boundary = shades of gray
            color = 255*(int)(distance);
        //}   
    }
   // if (nz < bailout) return 1; //still not escaping after iteration , rare
  // if (absdZ2 < absZ2) color= iExterior;  //includes escaping through 0 // som eiterior points are coded as exterior = error
   
  return color;
} //0.000007
 


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

int main(){

  unsigned int iX,iY, /* indices of 2D virtual array (image) = integer coordinate */
    i; /* index of 1D array  */
  double Cx, Cy;
 
  //int  ExtLastIteration; //
  unsigned char color;


  printf(" Setup \n");
  iXmax= iSide; /* height of image in pixels */
  iYmax  = iSide;
  iLength = iXmax*iYmax; /* number of pixels */
  CxMin = -2.25;
  CxMax = 0.75  ;
  CyMin = -dSide;
  CyMax = dSide;
  IterationMax = 1000 ;
 PixelWidth=  (CxMax-CxMin)/iXmax;
 PixelHeight = (CyMax-CyMin)/iYmax;
  ER2 = EscapeRadius*EscapeRadius;

  /* dynamic 1D array for colors ( shades of gray ) */
  unsigned char *data;
  data = malloc( iLength * sizeof(unsigned char) );
  if (data == NULL )
    {
      fprintf(stderr," Could not allocate memory");
      return 1;
    }
   
  printf(" compute color \n");
  for(iY=0;iY<iYmax;++iY){
    Cy=CyMin + iY*PixelHeight; /*  */
    printf("row %u from %u \n",iY, iYmax);   
    for(iX=0;iX<iXmax;++iX){
      Cx=CxMin + iX*PixelWidth;
      color=GiveColor(Cx, Cy, IterationMax, PixelWidth);
     
      i= f(iX,iY); /* compute index of 1D array from indices of 2D array */
      data[i]=color;  /* change the color */
     
    }
  }


 

   

  printf(" save  data array to the pgm file \n");
  unsigned int length = 30;
  char filename[length] ;
  snprintf(filename, length, "%.0fe%u%s", EscapeRadius,IterationMax, ".pgm");
  char *comment="#  ";/* comment should start with # */
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  /* save image to pgm file  */
  FILE * fp;     
  fp = fopen(filename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"P5\n %s\n %u\n %u\n %u\n",comment,iXmax,iYmax,MaxColorComponentValue);  /*write header to the file*/
  fwrite(data,iLength,1,fp);  /*write image data bytes to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);

  /* ---------- text file  -------------------------------------*/
  char tfilename[length] ;
  printf(" save  graphic data to the text file \n");
  snprintf(tfilename, length, "%.0fe%u%s",EscapeRadius, iXmax, ".txt");
   /* save info  to the file  */     
  fp = fopen(tfilename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"IterationMax = %d \n", IterationMax);
  fprintf(fp,"EscapeRadius = %f ER2 = %f \n", EscapeRadius, ER2);
 

  fprintf(fp,"\n" );  /* */
  fprintf(fp,"C plane : \n" );  /* */
  fprintf(fp,"dSide = %f \n", dSide);
  fprintf(fp,"PixelWidth = %f ; PixelHeight = %f \n", PixelHeight, PixelWidth);
  fprintf(fp," CxMin = %f \n", CxMin);  /* */
  fprintf(fp," CxMax = %f \n", CxMax);  /* */
  fprintf(fp," CyMin = %f \n", CyMin);  /* */
  fprintf(fp," CyMax = %f \n", CyMax);  /* */
  fprintf(fp,"\n" );
  //
  printf("File %s saved. \n", tfilename);
  fclose(fp);


  printf(" allways free memory \n ");
  free(data);
 
 

  return 0;
}
Logged
Adam Majewski
Fractal Lover
**
Posts: 221


WWW
« Reply #3 on: January 11, 2014, 09:36:51 AM »

Code:
distance = sqrt(absZ2/absdZ2)*log(absZ2); //
     distance = pow( 4.0*distance, 0.25 );
     if (distance<=1.0) color = (int)(255.0*distance);


Now it works.  Thx
[code
/*

  c console program, for CPU, one thread. numbers type : double
  It can be compiled and run under Linux, windows, Mac
  It needs gcc


draw :
* check 2 algorithms :
** binary escape time
** add boundary computed by DEM/M
* save it to the pgm file

 
 
 



  -----------------------------------------
  1.pgm file code is  based on the code of Claudio Rocchini
  http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
  create 8 bit color graphic file ,  portable gray map file = pgm
  see http://en.wikipedia.org/wiki/Portable_pixmap
  to see the file use external application ( graphic viewer)
  I think that creating graphic can't be simpler
  ---------------------------
  2. first it creates data array which is used to store color values of pixels,
  fills tha array with data and after that writes the data from array to pgm file.
  It alows free ( non sequential) acces to "pixels"
   
  -------------------------------------------
  Adam Majewski   fraktal.republika.pl
 
 

  to compile :
  gcc m.c  -lm -Wall -march=native
  to run ( Linux console) :
  time ./a.out




convert h6.650000m6650.pgm -resize 800x120 c.png
convert b6.650000m6650.pgm -resize 1600x240 cbig.png
 
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>



/* iXmax/iYmax = 1 */
unsigned int iSide = 1000; /* side of rectangle in pixels */
unsigned int iXmax; // ((int)(m*iSide)) /* height of image in pixels */
unsigned int iYmax ; //= iSide;
unsigned int iLength; // = (iXmax*iYmax) /* number of pixels */
/*  world ( double) coordinate */
double dSide =  1.5;
double CxMin ; // = 0.0;
double CxMax ; // =(m*dSide);
double CyMin ; //= dSide;
double CyMax ; // = dSide;
/* (CxMax-CxMin)/(CyMax-CyMin)==iXmax/iYmax = 1 */

unsigned int IterationMax; // = (iXmax*100) /* proportional to resolution of picture */

double PixelWidth; //=  ((CxMax-CxMin)/iXmax)
double PixelHeight ;//= ((CyMax-CyMin)/iYmax)

double CDistanceMax ; //= PixelWidth; /* proportional to pixel size */


/* fc(z) = z*z + c */



double EscapeRadius = 33.0;  /* radius of circle around origin; its complement is a target set for escaping points */
double ER2 ; //= (EscapeRadius*EscapeRadius)


/* colors = shades of gray from 0=black  to 255=white */
unsigned char iExterior = 245; /* exterior of Julia set */
unsigned char iBoundary = 0; /* border , boundary*/
unsigned char iInterior = 0;







unsigned int f(unsigned int _iX, unsigned int _iY)
/*
   gives position of point (iX,iY) in 1D array  ; uses also global variables
   it does not check if index is good  so memory error is possible
*/
{return (_iX + (iYmax-_iY-1)*iXmax );}








/*
 estimates distance from point c to nearest point in Julia  set
 for Fc(z)= z*z + c
 z(n+1) = Fc(zn) 
 this function is based on function  mndlbrot::dist  from  mndlbrot.cpp
 from program mandel by Wolf Jung (GNU GPL )
 http://www.mndynamics.com/indexp.html

Hyunsuk Kim  :
For Julia sets, z is the variable and c is a constant. Therefore df[n+1](z)/dz = 2*f[n]*f'[n] -- you don't add 1.

For the Mandelbrot set on the parameter plane, you start at z=0 and c becomes the variable. df[n+1](c)/dc = 2*f[n]*f'[n] + 1.
http://iquilezles.org/www/articles/distancefractals/distancefractals.htm

 */

// boolean Escape time and DEM/M in one loop 
unsigned char GiveColor( double Cx, double Cy, int iMax, double DistanceMax)
{

  // C = Cx + Cy* I = point of parameter c-plane

  int i; /* iteration number */

  double Zx, Zy; /* Z = Zx + Zy*I  point of dynamical plane */
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  double temp;
  double absZ2;

  // http://en.wikipedia.org/wiki/Complex_quadratic_polynomial
  // first derivative  of fc(zcr) with respect to c =  dZ = dfc(zcr)/dc = 2*Z*dZ  = dZx + dZy*I
  double dZx = 0.0;
  double dZy = 0.0;
  double absdZ2; // = abs(dZ)* abs(dZ) = dZx*dZx + dZy*dZy

  double distance;
  unsigned char color = iInterior;

 

  /* initial value of orbit  = critical point zcr = 0 */
  Zx=0.0;
  Zy=0.0;
  //
  Zx2 = Zx*Zx;
  Zy2 = Zy*Zy;
  absZ2= Zx*Zx + Zy*Zy;
  absdZ2= dZx*dZx + dZy*dZy;
  // iteration of critical point z= 0 on the dynamical z-plane       
  for (i=0; i<iMax; i++)
    { // check if not escaping : abs(z)>ER
      if (absZ2  > ER2 ) { color=  iExterior; break ; } // exterior when escapes
      //if (absdZ2 > 1e60) { return iInterior; } //  interior when derivative explodes
     

       // in the loop, the derivative should be calculated before the new z
      /* first derivative   zp = 2*z*zp  = xp + yp*i; */
      temp = 2*(Zx*dZx - Zy*dZy) + 1.0 ; /*  */
      dZy = 2*(Zx*dZy + Zy*dZx);
      dZx = temp;

      // z = fc(z) = z*z + c
      Zy=2*Zx*Zy + Cy;
      Zx=Zx2-Zy2 +Cx;
     
     

      // abs
      Zx2 = Zx*Zx;     
      Zy2 = Zy*Zy;
      absZ2= Zx2 + Zy2; // nz = x*x + y*y;
      absdZ2= dZx*dZx + dZy*dZy; //
    };

  if (i<iMax)
   {
     
     distance = sqrt(absZ2/absdZ2)*log(absZ2); //
     distance = pow( 4.0*distance, 0.25 );
     //if (distance > 0.0 && distance<DistanceMax*)
       //   { //color = ((int)(255.0*distance)) % 255;

            // d = sqrt( dot(z,z)/dot(dz,dz) )*log(dot(z,z));
        // do some soft coloring based on distance
    //d = pow( 4.0*d/zoom, 0.25 );
            //if (distance < (DistanceMax)) color = iBoundary ; // bondary = black
            // near boundary = shades of gray
            if (distance<=1.0) color = (int)(255.0*distance);
        //}   
    }
   // if (nz < bailout) return 1; //still not escaping after iteration , rare
  // if (absdZ2 < absZ2) color= iExterior;  //includes escaping through 0 // som eiterior points are coded as exterior = error
   
  return color;
} //0.000007
 


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

int main(){

  unsigned int iX,iY, /* indices of 2D virtual array (image) = integer coordinate */
    i; /* index of 1D array  */
  double Cx, Cy;
 
  //int  ExtLastIteration; //
  unsigned char color;


  printf(" Setup \n");
  iXmax= iSide; /* height of image in pixels */
  iYmax  = iSide;
  iLength = iXmax*iYmax; /* number of pixels */
  CxMin = -2.25;
  CxMax = 0.75  ;
  CyMin = -dSide;
  CyMax = dSide;
  IterationMax = 1000 ;
 PixelWidth=  (CxMax-CxMin)/iXmax;
 PixelHeight = (CyMax-CyMin)/iYmax;
  ER2 = EscapeRadius*EscapeRadius;

  /* dynamic 1D array for colors ( shades of gray ) */
  unsigned char *data;
  data = malloc( iLength * sizeof(unsigned char) );
  if (data == NULL )
    {
      fprintf(stderr," Could not allocate memory");
      return 1;
    }
   
  printf(" compute color \n");
  for(iY=0;iY<iYmax;++iY){
    Cy=CyMin + iY*PixelHeight; /*  */
    printf("row %u from %u \n",iY, iYmax);   
    for(iX=0;iX<iXmax;++iX){
      Cx=CxMin + iX*PixelWidth;
      color=GiveColor(Cx, Cy, IterationMax, PixelWidth);
     
      i= f(iX,iY); /* compute index of 1D array from indices of 2D array */
      data[i]=color;  /* change the color */
     
    }
  }


 

   

  printf(" save  data array to the pgm file \n");
  unsigned int length = 30;
  char filename[length] ;
  snprintf(filename, length, "%.0fe%u%s", EscapeRadius,IterationMax, ".pgm");
  char *comment="#  ";/* comment should start with # */
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  FILE * fp;     
  fp = fopen(filename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"P5\n %s\n %u\n %u\n %u\n",comment,iXmax,iYmax,MaxColorComponentValue);  /*write header to the file*/
  fwrite(data,iLength,1,fp);  /*write image data bytes to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);

  printf(" save  graphic data to the text file \n");
  char tfilename[length] ;
  snprintf(tfilename, length, "%.0fe%u%s",EscapeRadius, iXmax, ".txt");
  fp = fopen(tfilename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"IterationMax = %d \n", IterationMax);
  fprintf(fp,"EscapeRadius = %f ER2 = %f \n", EscapeRadius, ER2);
  fprintf(fp,"\n" );  /* */
  fprintf(fp,"C plane : \n" );  /* */
  fprintf(fp,"Sides  = %f x %f \n",CxMax- CxMin, CyMax - CyMin );
  fprintf(fp,"PixelWidth = %f ; PixelHeight = %f \n", PixelHeight, PixelWidth);
  fprintf(fp," CxMin = %f \n", CxMin);  /* */
  fprintf(fp," CxMax = %f \n", CxMax);  /* */
  fprintf(fp," CyMin = %f \n", CyMin);  /* */
  fprintf(fp," CyMax = %f \n", CyMax);  /* */
  fprintf(fp," center C = %f ; %f \n",CxMax -(CxMax-CxMin)/2.0, CyMax - (CyMax-CyMin)/2.0);  /* */
  fprintf(fp,"\n" );
  printf("File %s saved. \n", tfilename);
  fclose(fp);


  printf(" allways free memory \n ");
  free(data);
 
 

  return 0;
}
]
Logged
Pages: [1]   Go Down
  Print  
 
Jump to:  

Related Topics
Subject Started by Replies Views Last post
SHADERTOY not working Help & Support hgjf2 1 437 Last post September 21, 2013, 01:41:27 PM
by Nahee_Enterprises
GLSL error on shadertoy web sandbox Programming Mrz00m 0 2157 Last post November 03, 2013, 10:22:21 AM
by Mrz00m
How do i put colors on this 3d shadertoy code? 3D Fractal Generation Mrz00m 1 3127 Last post September 04, 2014, 01:08:29 PM
by cKleinhuis
shadertoy 3d peano? 3D Fractal Generation DarkBeam 6 4801 Last post March 08, 2015, 04:31:44 PM
by DarkBeam
shadertoy remixing Movies Showcase (Rate My Movie) cKleinhuis 6 6548 Last post August 02, 2015, 09:47:26 AM
by hgjf2

Powered by MySQL Powered by PHP Powered by SMF 1.1.21 | SMF © 2015, Simple Machines

Valid XHTML 1.0! Valid CSS! Dilber MC Theme by HarzeM
Page created in 0.18 seconds with 24 queries. (Pretty URLs adds 0.008s, 2q)