Fractals/Iterations in the complex plane/demm

This algorithm has 2 versions : exterior interior Compare it with version for dynamic plane and Julia set : DEM/J Exterior DEM/M right|thumb|Exterior DEM/M right|thumb|simple boundary with DEM/M right|thumb|Boundary with DEM/M and Sobel filter Distance Estimation Method for Mandelbrots set ( DEM/M ) estimates distance from point ( in the exterior of Mandelbrot set ) to nearest point in Mandelbrot set.A Cheritat wiki-draw: Mandelbrot_setfractalforums : m-set-distance-estimation Names the (directional) distance estimate formula It can be used to create B&W images from BOF :fractalforums discussion : How are B&W images from "Beauty of Fractals" created? map 41 on page 84 map 43 on page 85 an unnumbered plot on page 188 types analytic DE = DEM = true DE Claude DE = fake DE ( FDE): FDE=1/(g*log(2)) when g>0 and undefined when g=0 (i.e., in interior)fractalforums.org : m-brot-distance-estimation-versus-claudes-fake-de Math formule Math formula :Milnor (in Dynamics in One Complex Variable, appendix G)Heinz-Otto Peitgen (Editor, Contributor), Dietmar Saupe (Editor, Contributor), Yuval Fisher (Contributor), Michael McGuire (Contributor), Richard F. Voss (Contributor), Michael F. Barnsley (Contributor), Robert L. Devaney (Contributor), Benoit B. Mandelbrot (Contributor) : The Science of Fractal Images. Springer; 1 edition (July 19, 1988), page 199 distance rendering for fractals by ińigo quilez where : is a first derivative of with respect to c which can be found by iteration starting with and then replacing at every consecutive step Escape radius "you need to change "If mag > 2" to increase the escape radius. Attached is a comparison, showing that the "straps" reduce when the radius is increased, no harm going huge like R = 1e6 or so because once it escapes 2 it grows very fast. The bigger the radius, the more accurate the distance estimate." Claude fractalforums : mandelbrot-distance-estimation-problem Basically in testing magnitude squared against bailout if your bailout is 10^4 then your DE estimates are accurate to 2 decimal figures i.e. if your DE is 0.01 then it's only accurate to +/-0.0001, if your bailout is 10^6 then your DE estimates are accurate to 3 decimal figures (David Makin in fractalforums.com) "For exterior distance estimation, you need a large escape radius, eg 100100. Iteration count limit is arbitrary, with a finite limit some pixels will always be classified as "unknown"." Claude Heiland-Allenfractalforums.org : list-of-numbers-with-fixed-digit-length-in-the-mandelbrot-set Comparison The Beauty of Fractals The Science of Fractal Images, page 198, Distance Estimator by Robert P. MunafoDistance Estimator by Robert P. Munafo Pseudocode by Claude Heiland-Allenmath.stackexchange question: how-to-draw-a-mandelbrot-set-with-the-connecting-filaments-visible foreach pixel c while not escaped and iteration limit not reached dc := 2 * z * dc + 1 z := z^2 + c de := 2 * |z| * log(|z|) / |dz| d := de / pixel_spacing plot_grey(tanh(d)) Code Mandelbrot book double _Complex m_exterior_distance(int N, double R, double _Complex c) { double _Complex z = 0; double _Complex dc = 0; for (int n = 0; n < N; ++n) { if (cabs(z) > R) return 2 * z * log(cabs(z)) / dc; dc = 2 * z * dc + 1; z = z * z + c; } return -1; } Analysis of code from BOF "The Beauty of Fractals" gives an almost correct computer program for the distance estimation shown in the right image. A possible reason that that method did not gain ground is that the procedure in this program is seriously flawed: The calculation of is performed (and completed) before the calculation of , and not after as it ought to be ( uses , not ). For the successive calculation of , we must know (which in this case is 2). In order to avoid the calculation of (k = 0, 1, 2, ...) again, this sequence is saved in an array. Using this array, is calculated up to the last iteration number, and it is stated that overflow can occur. If overflow occurs the point is regarded as belonging to the boundary (the bail-out condition). If overflow does not occur, the calculation of the distance can be performed. Apart from it being untrue that overflow can occur, the method makes use of an unnecessary storing and repetition of the iteration, making it unnecessarily slower and less attractive. The following remark in the book is nor inviting either: "It turns out that the images depend very sensitively on the various choices" (bail-out radius, maximum iteration number, overflow, thickness of the boundary and blow-up factor). Is it this nonsense that has got people to lose all desire for using and generalizing the method? " Gertbuschmann So : in the loop, derivative should be calculated before the new z "I'm finally generating DEM (Distance Estimate Method) based images that I'm happy with. It turns out that my code had a couple of bugs in it. The new code runs reasonably fast, even with all the extra math to compute distance estimates. The algorithm in S of F I ("The Science of Fractal Images" uses a sharp cut-off from white to black when pixels get close enough to the set, but I find this makes for jagged looking plots. I've implemented a non-linear color gradient that works pretty well. For pixels where their DE value is threshold>=DE>=0, I scale the value of DE to 1>=DE>=0, then do the calculation: gradient_value = 1 - (1-DE)^n, ("^n" means to the n power. I wish I could use superscripts!) "n" is an adjustment factor that lets me change the shape of the gradient curve. The value "gradient_value" determines the color of the pixel. If it's 0, the pixel is colored at the starting color (white, for a B&W plot.) If "gradient_value" is 1, the pixel gets the end color (e.g. black) For values of n>1, this gives a rapid change in color for pixels that are far from the set, and a slower change in value as DE approaches 0. For values of n<1, it gives a slow color change for pixels far from the set, and rapid change close to the set. For n=1, I get a linear color gradient. The non-linear gradient lets me use the DE value to anti-alias my plots. By adjusting my threshold value and my adjustment value, I can get good looking results for a variety of images." Duncan C Fractal forum : How are B&W images from "Beauty of Fractals" created? "All our DE formulas above are only approximations – valid in the limit n→∞, and some also only for point close to the fractal boundary. This becomes very apparent when you start rendering these structures – you will often encounter noise and artifacts. Multiplying the DE estimates by a number smaller than 1, may be used to reduce noise (this is the Fudge Factor in Fragmentarium). Another common approach is to oversample – or render images at large sizes and downscale." Mikael Hvidtfeldt Christensen Distance Estimated 3D Fractals (V): The Mandelbulb & Different DE Approximations by Mikael Hvidtfeldt Christensen Example code ultrafractal Two algorithm in two loops Distance estimation to the Mandelbrot set by Inigo Quilez in shadertoy : Distance estimation to the Mandelbrot set by Inigo Quilez Here is Java function. It computes in one loop both : iteration and derivative . If (on dynamic plane) critical point : does not escapes to infinity ( bailouts), then it belongs to interior of Mandelbrot set and has colour maxColor, escapes then it maybe in exterior or near boundary. Its colour is "proportional" to distance between the point c and the nearest point in the Mandelbrot set. It uses also colour cycling ( (int)R % maxColor ). // Java function by Evgeny Demidov from http://www.ibiblio.org/e-notes/MSet/DEstim.htm // based on code by Robert P. Munafo from http://www.mrob.com/pub/muency/distanceestimator.html public int iterate(double cr, double ci, double K, double f) { double Cr,Ci, I=0, R=0, I2=I*I, R2=R*R, Dr=0,Di=0,D; int n=0; if (f == 0){ Cr = cr; Ci = ci;} else{ Cr = ci; Ci = cr;} do { D = 2*(R*Dr - I*Di) +1; Di = 2*(R*Di + I*Dr); Dr = D; I=(R+R)*I+Ci; R=R2-I2+Cr; R2=R*R; I2=I*I; n++; } while ((R2+I2 < 100.) && (n < MaxIt) ); if (n == MaxIt) return maxColor; // interior else{ // boundary and exterior R = -K*Math.log(Math.log(R2+I2)*Math.sqrt((R2+I2)/(Dr*Dr+Di*Di))); // compute distance if (R < 0) R=0; return (int)R % maxColor; }; } Here is cpp function. It gives integer index of color as an output. /* this function is from program mandel ver 5.10 by Wolf Jung see file mndlbrot.cpp http://www.mndynamics.com/indexp.html It is checked first (in pixcolor) that the point escapes before calling this function. So we do not compute the derivative and the logarithm when it is not escaping. This would not be a good idea if most points escape anyway. */ int mndlbrot::dist(double a, double b, double x, double y) { uint j; double xp = 1, yp = 0; // zp = xp+yp*i = 1 ; derivative double nz, nzp; // iteration for (j = 1; j <= maxiter; j++) { // zp nz = 2*(x*xp - y*yp); yp = 2*(x*yp + y*xp); xp = nz; //zp = 2*z*zp; on the dynamic plane // if sign is positive it is parameter plane, if negative it is dynamic plane. if (sign > 0) xp++; //zp = 2*z*zp + 1 ; on the parameter plane //z = z*z + c; nz = x*x - y*y + a; y = 2*x*y + b; x = nz; // nz = x*x + y*y; nzp = xp*xp + yp*yp; // 2 conditions for stopping the iterations if (nzp > 1e40 || nz > bailout) break; } // 5 types of points but 3 colors /* not escaping, rare Is it not simply interior ??? This should not be necessary. I do not know why I included it, because in this case pixcolor should not call dist. If you do not have pixcolor before, you should return 0 here. */ if (nz < bailout) return 1; // not escaping, rare /* If The Julia set is disconnected and the orbit of z comes close to z = 0 before escaping, nzp will be small */ if (nzp < nz) return 10; // includes escaping through 0 // compute estimated distance = x x = log(nz); x = x*x*nz / (temp[1]*nzp); //4*square of dist/pixelwidth if (x < 0.04) return 1; // exterior but very close to the boundary if (x < 0.24) return 9; // exterior but close to the boundary return 10; //exterior : escaping and not close to the boundary } //dist Two algorithms in one loop distance rendering for fractals by ińigo quilez in c++distance rendering for fractals by ińigo quilez C source code for CPU one thread DEM/M 8 bit - Click on the right to view/* 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 the array with data and after that writes the data from array to pgm file. It allows free ( non sequential) access 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 = 255; /* 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 double EscapeRadius = 33.0; */ // boolean Escape time and DEM/M in one loop double GiveDistance( 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; /* 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 ) break ; // exterior when escapes //if (absdZ2 > 1e60) { i=iMax; } // 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; // }; // compute distance if (i<iMax) // exterior { // based on https://www.shadertoy.com/view/lsX3W4 shadertoy : Distance estimation to the Mandelbrot set by Inigo Quilez distance = sqrt(absZ2/absdZ2)*log(absZ2); // distance = pow( 4.0*distance, 0.25 ); } else distance=0.0; // interior // 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 distance; } //0.000007 // color is proportional to distance between point c and nearest point in Mandelbrot set unsigned char GiveColor( double Cx, double Cy, int iMax, double DistanceMax) { double distance ; unsigned char color ; distance = GiveDistance( Cx, Cy, iMax, DistanceMax); if (distance > 0.0) { if (distance<DistanceMax*333.0) // = clamp(distance, 0.0, 1.0) to remove level sets effect !!!!! color = (int)(255.0*distance); // boundary and near boundary = shades of gray else color = iExterior; // exterior far away from boundary } else color = iInterior; return color; } /* --------------------------------------------------------------------------------------------------------- */ 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 for every pixel : scan c-plane \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,iXmax , ".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,"dWidth = %f ; dHeight = %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 of image : 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; } C source code for CPU one thread DEM/M 32 bitcolor - Click on the right to view/* 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 ppm file ----------------------------------------- 1.pgm file code is based on the code of Claudio Rocchini http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg create 24 bit color graphic file , portable pixel map file = ppm 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 the array with data and after that writes the data from array to pgm file. It allows free ( non sequential) access 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> double g=100.0; /* 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 = 255; /* 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)*3 );} // ((iYmax-iY-1)*iXmax+iX)*3; /* based on Delphi function by Witold J.Janik */ void GiveRainbowColor(double position,unsigned char Color32[]) { /* if position > 1 then we have repetition of colors it maybe useful */ if (position>1.0){if (position-(int)position==0.0)position=1.0; else position=position-(int)position;} unsigned char nmax=6; /* number of color segments */ double m=nmax* position; int n=(int)m; // integer of m double f=m-n; // fraction of m unsigned char t=(int)(f*255); switch( n){ case 0: { Color32[0] = 255; Color32[1] = t; Color32[2] = 0; break; }; case 1: { Color32[0] = 255 - t; Color32[1] = 255; Color32[2] = 0; break; }; case 2: { Color32[0] = 0; Color32[1] = 255; Color32[2] = t; break; }; case 3: { Color32[0] = 0; Color32[1] = 255 - t; Color32[2] = 255; break; }; case 4: { Color32[0] = t; Color32[1] = 0; Color32[2] = 255; break; }; case 5: { Color32[0] = 255; Color32[1] = 0; Color32[2] = 255 - t; break; }; default: { Color32[0] = 255; Color32[1] = 0; Color32[2] = 0; break; }; }; // case } /* 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 double EscapeRadius = 33.0; */ // boolean Escape time and DEM/M in one loop double GiveDistance( 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; /* 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 ) break ; // exterior when escapes //if (absdZ2 > 1e60) { i=iMax; } // 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; // }; // compute distance if (i<iMax) // exterior { distance = sqrt(absZ2/absdZ2)*log(absZ2); // distance = pow( g*distance, 0.25 ); } else distance=0.0; // interior // 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 distance; } //0.000007 // gives width in estimated distance of near boundary strip around Mandelbrot set // true boundary is = DistanceMax*0.1 double GiveDistanceMax(double PixelWidth, int iSide) { return PixelWidth*iSide*0.0333; // = 0.1 } // color is proportional to distance between point c and nearest point in Mandelbrot set int ComputeColor32( double Cx, double Cy, int iMax, double DistanceMax, unsigned char color32[3] ) { double distance ; distance = GiveDistance( Cx, Cy, iMax, DistanceMax); if (distance > 0.0) { if (distance<DistanceMax) // boundary { // boundary color32[0] = 255-(int)(255.0*distance); color32[1] = 255-(int)(255.0*distance); color32[2] = 255-(int)(255.0*distance); } else { // exterior far away from boundary GiveRainbowColor(distance, color32); // gradient of 32 rgb color } } else { // interior color32[0] = 0; color32[1] = 0; color32[2] = 0; } return 0; } /* --------------------------------------------------------------------------------------------------------- */ 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 color32[3]; printf(" Setup \n"); iXmax = iSide; /* height of image in pixels */ iYmax = iSide; iLength = iXmax*iYmax*3; /* number of pixels* number of bytes per color */ CxMin = -2.25; CxMax = 0.75 ; CyMin = -dSide; CyMax = dSide; IterationMax = 1000 ; PixelWidth= (CxMax-CxMin)/iXmax; PixelHeight = (CyMax-CyMin)/iYmax; ER2 = EscapeRadius*EscapeRadius; CDistanceMax = GiveDistanceMax( PixelWidth, iSide); /* 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 for every pixel : scan c-plane \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; ComputeColor32(Cx, Cy, IterationMax, CDistanceMax, color32); i= f(iX,iY); /* compute index of 1D array from indices of 2D array */ /* change the color */ data[i] = color32[0]; data[i+1] = color32[1]; data[i+2] = color32[2]; } } printf(" save data array to the pgm file \n"); unsigned int length = 30; char filename[length] ; snprintf(filename, length, "%fg%.0fe%u%s",g, EscapeRadius,iXmax , ".ppm"); 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,"P6\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, "%fg%.0fe%u%s",g, 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,"dWidth = %f ; dHeight = %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 of image : C = %f ; %f \n",CxMax -(CxMax-CxMin)/2.0, CyMax - (CyMax-CyMin)/2.0); /* */ fprintf(fp,"\n" ); fprintf(fp,"g= %f used for distance = pow( g*distance, 0.25 ); \n", g); printf("File %s saved. \n", tfilename); fclose(fp); printf(" always free memory \n "); free(data); return 0; } Function GiveDistance(xy2,eDx,eDy:extended):extended; begin result:=2*log2(sqrt(xy2))*sqrt(xy2)/sqrt(sqr(eDx)+sqr(eDy)); end; //------------------------------------ Function PointIsInCardioid (Cx,Cy:extended):boolean; //Hugh Allen // http://homepages.borland.com/ccalvert/Contest/MarchContest/Fractal/Source/HughAllen.zip var DeltaX,DeltaY:extended; // PDeltyX,PDeltyY:extended; // ZFixedX,ZFixedY:extended; begin result:=false; // cardioid checkig - thx to Hugh Allen //sprawdzamy Czy punkt C jest w głównej kardioidzie //Cardioid in squere :[-0.75,0.4] x [ -0.65,0.65] if InRange(Cx,-0.75,0.4)and InRange(Cy,-0.65,065) then begin // M1= all C for which Fc(z) has attractive( = stable) fixed point // znajdyjemy punkt staly z: z=z*z+c // czyli rozwiazujemy rownanie kwadratowe // zmiennej zespolonej o wspolczynnikach zespolonych // Z*Z - Z + C = 0 //Delta:=1-4*a*c; Delta i C sa liczbami zespolonymi DeltaX:=1-4*Cx; DeltaY:=-4*Cy; // Pierwiastek zespolony z delty CmplxSqrt(DeltaX,DeltaY,PDeltyX,PDeltyY); // obliczmy punkt staly jeden z dwóch, ten jest prawdopodobnie przycišgajšcy ZFixedX:=1.0-PDeltyX; //0.5-0.5*PDeltyX; ZFixedY:=PDeltyY; //-0.5*PDeltyY; // jesli punkt stały jest przycišgajšcy // to należy do M1 If (ZfixedX*ZFixedX + ZFixedY*ZFixedY)<1.0 then result:=true; // ominięcie iteracji M1 przyspiesza 3500 do 1500 msek end; // if InRange(Cx ... end; //------------------------------------ Function PointIsInComponent (Cx,Cy:extended):boolean; //Hugh Allen // http://homepages.borland.com/ccalvert/Contest/MarchContest/Fractal/Source/HughAllen.zip var Dx:extended; begin result:=false; // czy punkt C nalezy do koła na lewo od kardioidy // circle: center = -1.0 and radius 1/4 dx:=Cx+1.0; if (Dx*Dx+Cy*Cy) < 0.0625 then result:=true; end; //------------------------------ Procedure DrawDEM_DazibaoTrueColor; // draws Mandelbrot set in black and its complement in true color // see http://ibiblio.org/e-notes/MSet/DEstim.htm // by Evgeny Demidov // // see also //http://www.mandelbrot-dazibao.com/Mset/Mset.htm // translation ( with modification) of Q-Basic program: // http://www.mandelbrot-dazibao.com/Mset/Mdb3.bas // // see also my page http://republika.pl/fraktal/mset_dem.html var iter:integer; iY,iX:integer; eCy ,eCx:extended; // C:=eCx + eCy*i eX,eY:extended; // Zn:=eX+eY*i eTempX,eTempY:extended; eX2,eY2:extended; //x2:=eX*eX; y2:=eY*eY; eXY2:extended; // xy2:=x2+y2; eXY4:extended; eTempDx,eTempDy:extended; eDx,eDy:extended; // derivative distance:extended; color:TColor; begin //compute bitmap for iY:= iYmin to iYMax do begin eCy:=Convert_iY_to_eY(iY); for iX:= iXmin to iXmax do begin eCx:=Convert_iX_to_eX(iX); If not PointIsInCardioid (eCx,eCy) and not PointIsInComponent(eCx,eCy) then begin // Z0:=0+0*i eX:=0; eY:=0; eTempX:=0; eTempY:=0; // eX2:=0; eY2:=0; eXY2:=0; // eDx:=0; eDy:=0; eTempDx:=0; eTempDy:=0; // iter:=0; // iteration of Z ; Z= Z*z +c while ((iter<IterationMax) and (eXY2<=BailOut2)) do begin inc(iter); // eTempY:=2*eX*eY + eCy; eTempX:=eX2-eY2 + eCx; // eX2:=eTempX*eTempX; eY2:=eTempY*eTempY; // eTempDx:=1+2*(eX*eDx-eY*eDy); eTempDy:=2*(eX*eDY+eY*eDx); // eXY2:=eX2+eY2; // eX:=ETempX; eY:=eTempY; // eDx:=eTempDx; eDy:=eTempDy; end; // while // drawing procedure if (iter<IterationMax) then begin distance:= GiveDistance(eXY2,eDx,eDy); color:=Rainbow(1,500,Abs(Round(100*Log10(distance)) mod 500)); with Bitmap1.FirstLine[iY*Bitmap1.LineLength+iX] do begin B := GetBValue(color); G := GetGValue(color); R := GetRValue(color); //A := 0; end; // with FirstLine[Y*LineLength+X] end // if (iter<IterationMax) then else with Bitmap1.FirstLine[iY*Bitmap1.LineLength+iX] do begin B := 0; G := 0; R := 0; //A := 0; end; //--- end of drawing procedure --- end // If not PointIsInCardioid ... then else with Bitmap1.FirstLine[iY*Bitmap1.LineLength+iX] do begin B := 0; G := 0; R := 0; //A := 0; end; //--- If not PointIsInCardioid ... end; // for iX end; // for iY end; //------------------------------------------------------ Modifications / optimisations Creating DEM images can be improved by use of : the Koebe 1/4-theoremdistance estimation by Claude Heiland-Allen gray scale based on distance distance rendering for fractals by ińigo quilezfractalforums gallery by Pauldelbrot Fractal forum : Relationship between bailout and accuracy of analytical DE weight in DEM/M if you see double contur of main antenna : increasing resolution does't solve the problem add one pixel to height or change CyMax or CyMin ( add one pixel size). equalisation Distance Estimate Equalisationfractalforums.org : histogram-de-coloring Program exrdeeqDistance-Estimate-Equalisation by Claude Heiland-Allen reads DEX and DEY channels outputs Y channel DE less than 0 becomes 0 DE greater than 1 becomes 1 otherwise Y is the index of DE in the sorted array, scaled to the range 0 to 1 exrdeeq in.exr out.exr Examples B of F map 43 DEM by Duncan ChampneyB of F map 43 DEM by Duncan Champney This plot is centered on -0.7436636773584340, 0.1318632143960234i at a width of about 6.25E-11 Smily Kaleidoscope by Duncan Champney Smily Kaleidoscope by Duncan Champney height=800 width=800 max_iterations=10000 center_r=-1.74920463346 center_i=-0.000286846601561 r_width=3.19E-10 Interior distance estimation Estimates distance from limitly periodic point ( in the interior of Mandelbrot set ) to the boundary of Mandelbrot set. Descriptions :Interior and exterior distance bounds for the Mandelbrot by Albert Lobo Cusidó Book : The Science of Fractal Images, Springer Verlag, Tokyo, Springer, New York, 1988 by Heinz-Otto Peitgen , D. Saupe, page 294 Interior and exterior distance bounds for the Mandelbrot by Albert Lobo Cusidó interior distance estimate by R Munafo mandelbrot-numerics library : m_d_interior by Claude thumb|Pixels colored according to the estimated interior distance ‎// https://mathr.co.uk/mandelbrot/book-draft/#interior-distance // Claude Heiland-Allen #include <complex.h> #include <math.h> double cnorm(double _Complex z) { return creal(z) * creal(z) + cimag(z) * cimag(z); } double m_interior_distance (double _Complex z0, double _Complex c, int p) { double _Complex z = z0; double _Complex dz = 1; double _Complex dzdz = 0; double _Complex dc = 0; double _Complex dcdz = 0; for (int m = 0; m < p; ++m) { dcdz = 2 * (z * dcdz + dz * dc); dc = 2 * z * dc + 1; dzdz = 2 * (dz * dz + z * dzdz); dz = 2 * z * dz; z = z * z + c; } return (1 - cnorm(dz)) / cabs(dcdz + dzdz * dc / (1 - dz)); } double m_distance(int N, double R, double _Complex c) { double _Complex dc = 0; double _Complex z = 0; double m = 1.0 / 0.0; int p = 0; for (int n = 1; n <= N; ++n) { dc = 2 * z * dc + 1; z = z * z + c; if (cabs(z) > R) return 2 * cabs(z) * log(cabs(z)) / cabs(dc); if (cabs(z) < m) { m = cabs(z); p = n; double _Complex z0 = m_attractor(z, c, p); double _Complex w = z0; double _Complex dw = 1; for (int k = 0; k < p; ++k) { dw = 2 * w * dw; w = w * w + c; } if (cabs(dw) <= 1) return m_interior_distance(z0, c, p); } } return 0; } Math description The estimate is given by : where is the period = length of the periodic orbit is the point from parameter plane to be estimated is the complex quadratic polynomial is the -fold iteration of is any of the points that make the periodic orbit ( limit cycle ) 4 derivatives of evaluated at : First partial derivative with respect to z It must be computed recursively by applying the chain rule : Starting points : First iteration : Second iteration : General rule for p iteration : First partial derivative with respect to c It must be computed recursively by applying the chain rule : Starting points : First iteration : Second iteration : General rule for p iteration : Second order partial derivative with respect to z It must be computed : together with : and recursively by applying the chain rule Starting points : First iteration : Second iteration : General rule for p iteration : Second order mixed partial derivative Algorithm Steps choose c check if critical point for given c is periodic ( interior ) or not ( exterior or near boundary or to big period ) if not periodic do not use this algorithm ( use exterior version) if periodic compute period and periodic point using compute distance using p iteration Computing period and periodic point Computing distance Code Java code by Albert Lobo CusidóInterior and exterior distance bounds for the Mandelbrot by Albert Lobo Cusidó Haskell code by ClaudeInterior coordinates in the Mandelbrot set by Claude and image Fractal forum : Coloring points inside the Mandelbrot set processing code MandelbrotDE in processing by ttoinou and descriptionfractalforums : classic-mandelbrot-with-distance-and-gradient-for-coloring by tit_toinou Problems There are two practical problems with the interior distance estimate: first, we need to find precisely, second, we need to find precisely. The problem with is that the convergence to by iterating requires, theoretically, an infinite number of operations. The problem with period is that, sometimes, due to rounding errors, a period is falsely identified to be an integer multiple of the real period (e.g., a period of 86 is detected, while the real period is only 43=86/2). In such case, the distance is overestimated, i.e., the reported radius could contain points outside the Mandelbrot set. For interior distance estimation, you need the period, then a number (maybe 10 or so should usually be enough) Newton steps to find the limit cycle. the interior checking code absolutely requires the reference to be at the nucleus of the island (not any of its child discs, and certainly not some random point) Optimisation Analogous to the exterior case, once b is found, we know that all points within the distance of b/4 from c are inside the Mandelbrot set. Adaptive super-sampling using distance estimate Adaptive super-sampling using distance estimate by Claude Heiland-Allen References Distance Estimation for Hybrid Escape Time Fractals by Claude Heiland-Allen 2020-04-23 Milnor "Self-similarity and hairiness in the Mandelbrot set" in Computers in Geometry and Topology, 1989.