File:Parabolic Julia set for internal angle 1 over 5 - new method.png
Size of this preview: 600 × 600 pixels. Other resolutions: 240 × 240 pixels | 480 × 480 pixels | 768 × 768 pixels | 1,024 × 1,024 pixels | 2,001 × 2,001 pixels.
Original file (2,001 × 2,001 pixels, file size: 33 KB, MIME type: image/png)
 |
This is a file from the Wikimedia Commons |
Summary
| Description |
English: Parabolic Julia set for fc(z) = z^2 + c where c is on boundary of main cardioid with internal angle 1 over 5. Drawing : new method |
| Date | |
| Source | Own work |
| Author | Adam majewski |
Algorithm
First step is to divide points z of complex dynamical plane into 2 subsets ( using bailout test; see function GiveColor ) :
- escaping points = exterior
- non-escaping points interior and boundary = filled Julia set
iter =GiveLastIteration(Zx, Zy );
if (iter< iterMax ) // point escapes
{ color = 245; } // exterior
else // Filled Julia Set
Second step is to divide points of interior into 5 subsets using :
- target set
- fall-in test
Target set features :
- center = alfa.
- radius = AR
- is divided into 5 subsets ( petals) by 5-arm star
- all points of interior fall into fixed point alfa
Fall-in test : iterate z and check when point z is inside target set.
Note that test is done after every not after every
dX=Zx-_ZAx;
dY=Zy-_ZAy;
d=dX*dX+dY*dY; // d = square of distance from zn to alfa
while ( d>_AR2) // check if z is inside target
{
for(i=0;i<period ;++i) // iMax = period !!!!
Find in which subset of target set point z falls using relative angle ( see function GiveNumber ) :
if ((0.0<angle) && (angle<0.2)) return 0;
if ((0.2<angle) && (angle<0.4)) return 1;
if ((0.4<angle) && (angle<0.6)) return 2;
if ((0.6<angle) && (angle<0.8)) return 3;
return 4; // angle > 0.8
The star is turned [1] so
angle =GiveTurn(dX +dY*I) - 0.17298227404701; // this turn offset is computed with Maxima CAS file
So range of angle in which point zn falls into target set gives a color of subset of interior of Julia set :
unsigned char colors[]={255,230,180, 160,140}; // number of colors >= period
// ....
int number =GiveNumber(Zx,Zy, Cx, Cy, period, AR2, creal(alfa), cimag(alfa));
color = colors[number];
Third step : when all points of plane are colored apply edge detection with Sobel filter ( see function AddBoundaries )
Forth step : save memory array to pgm file
Five step : convert pgm to png using Image Magic :
convert a.pgm a.png
Licensing
I, the copyright holder of this work, hereby publish it under the following licenses:
This file is licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license.
- You are free:
- to share – to copy, distribute and transmit the work
- to remix – to adapt the work
- Under the following conditions:
- attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
- share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.
|
Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of the license is included in the section entitled GNU Free Documentation License. |
You may select the license of your choice.
Maxima CAS src code
This program computes turn offset = angle in turns. This angle equal to rotation of 5-arm star with center in alfa fixed point.
It works here but for other values ( period , q) do not. Help is welcome to improve this code !!!!
kill(all); remvalue(all); /* ------- functions ------------ */ f(z,c):=z*z+c; fn(p, z, c) := if p=0 then z elseif p=1 then f(z,c) else f(fn(p-1, z, c),c); /* gives parameter c for complex quadratic polynomial */ GiveC(angle,radius):= ( [w,c], w:radius*%e^(%i*angle*2*%pi), /* point of circle */ /* conformal map from circle to cardioid ( boundary of period 1 component of Mandelbrot set */ c:w/2-w*w/4, float(rectform(c)) /* point on boundary of period 1 component of Mandelbrot set */ )$ /* converts angle in radians in range -Pi to Pi to turns */ GiveTurn(a):= ( [a,t], if a<0 then a:a+2*%pi, /* from range (-Pi,Pi) to range (0,2Pi) */ t:a/(2*%pi) /* from radians to turns */ )$ /* http://num.math.uni-goettingen.de/~summer/cheritat.pdf PARABOLIC IMPLOSION. A MINI-COURSE by ARNAUD CHERITAT Let Czk be the next term in the expansion of fq: fq(z) = z + Czk + : : : We defne attracting axes as the k-1 half lines whose union is the solution to (C*z^6)/z < 0 ; and the repelling axes are for (C*z^6)/z > 0 ; */ GiveTermC(p,q):= ([fq, c, C], fq:fn(q,z,c), C:coeff(rat(fq),z,q+1), c:GiveC(p/q,1), /* parabolic point : rational angle and radius=1.0 */ C:rectform(float(ev(C))) )$ GiveTurnOffset(p,q):= ( [C,a,t,offset], C:GiveTermC(p,q), a:carg(C), t:GiveTurn(a), tt:t/q, /* turn offset */ tt:float(rectform(tt)) )$ compile(all); /* --------- const ----------- */ p:1; q:5; /* m:exp(2*%i*%pi*p/q); multiplier */ /* --------------- main ----------------------------- */ TurnOffset : GiveTurnOffset(p,q); /* // this turn offset is computed with Maxima CAS filename // star is slightly turned so this offset // theory : Complex Dynamics by Lennart Carleson,Theodore Williams Gamelin page 40 // TurnOffset = carg(a)/q // where a = coeff( z^(q+1)) */
C src code
/*
Adam Majewski
fraktal.republika.pl
c console progam using
* symmetry
* openMP
draw julia sets
gcc t.c -lm -Wall -fopenmp -march=native
time ./a.out
iHeight =1000
iterMax = 1000
AR = PixelWidth*15.0; // if you will increase iHeight
then increase multiplier or the time of computations will be very long
real 3m31.345s
File new2001.pgm saved.
Cx = 0.356763
Cy = 0.328582
alfax = 0.154508
alfay = 0.475528
distorsion of image = 1.000000
real 407m35.605s
*/
#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; needs also -fopenmp
/* --------------------------------- global variables and consts ------------------------------------------------------------ */
// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1
unsigned int ix, iy; // var
unsigned int ixMin = 0; // Indexes of array starts from 0 not 1
unsigned int ixMax ; //
unsigned int iWidth ; // horizontal dimension of array
unsigned int ixAxisOfSymmetry ; //
unsigned int iyMin = 0; // Indexes of array starts from 0 not 1
unsigned int iyMax ; //
unsigned int iyAxisOfSymmetry ; //
unsigned int iyAbove ; // var, measured from 1 to (iyAboveAxisLength -1)
unsigned int iyAboveMin = 1 ; //
unsigned int iyAboveMax ; //
unsigned int iyAboveAxisLength ; //
unsigned int iyBelowAxisLength ; //
unsigned int iHeight = 2000; // odd number !!!!!! = (iyMax -iyMin + 1) = iyAboveAxisLength + iyBelowAxisLength +1
// The size of array has to be a positive constant integer
unsigned int iSize ; // = iWidth*iHeight;
// memmory 1D array
unsigned char *data;
unsigned char *edge;
// unsigned int i; // var = index of 1D array
unsigned int iMin = 0; // Indexes of array starts from 0 not 1
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
/* world ( double) coordinate = dynamic plane */
const double ZxMin=-1.5;
const double ZxMax=1.5;
const double ZyMin=-1.5;
const double ZyMax=1.5;
double PixelWidth; // =(ZxMax-ZxMin)/iXmax;
double PixelWidth2; // = PixelWidth*PixelWidth;
double PixelHeight; // =(ZyMax-ZyMin)/iYmax;
double distanceMax;
double ratio ;
double TwoPi=2*M_PI;
// complex numbers of parametr plane
double Cx; // c =Cx +Cy * i
double Cy;
double complex c; //
double complex alfa; // alfa fixed point alfa=f(alfa)
unsigned long int iterMax = 1000; //iHeight*100;
double ER = 2.0; // Escape Radius for bailout test
double ER2;
double AR ; // radius around attractor
double AR2; // =AR*AR;
unsigned int period = 5; // period of secondary component joined by root point
unsigned int denominator;
double InternalAngle;
unsigned char colors[]={255,230,180, 160,140}; // number of colors >= period
/* ------------------------------------------ functions -------------------------------------------------------------*/
/* find c in component of Mandelbrot set
uses code by Wolf Jung from program Mandel
see function mndlbrot::bifurcate from mandelbrot.cpp
http://www.mndynamics.com/indexp.html
*/
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int period)
{
//0 <= InternalRay<= 1
//0 <= InternalAngleInTurns <=1
double t = InternalAngleInTurns *2*M_PI; // from turns to radians
double R2 = InternalRadius * InternalRadius;
//double Cx, Cy; /* C = Cx+Cy*i */
switch ( period ) // of component
{
case 1: // main cardioid
Cx = (cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4;
Cy = (sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4;
break;
case 2: // only one component
Cx = InternalRadius * 0.25*cos(t) - 1.0;
Cy = InternalRadius * 0.25*sin(t);
break;
// for each period there are 2^(period-1) roots.
default: // higher periods : to do
Cx = 0.0;
Cy = 0.0;
break; }
return Cx + Cy*I;
}
/*
http://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings
z^2 + c = z
z^2 - z + c = 0
ax^2 +bx + c =0 // ge3neral for of quadratic equation
so :
a=1
b =-1
c = c
so :
The discriminant is the d=b^2- 4ac
d=1-4c = dx+dy*i
r(d)=sqrt(dx^2 + dy^2)
sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx +- sy*i
x1=(1+sqrt(d))/2 = beta = (1+sx+sy*i)/2
x2=(1-sqrt(d))/2 = alfa = (1-sx -sy*i)/2
alfa : attracting when c is in main cardioid of Mandelbrot set, then it is in interior of Filled-in Julia set,
it means belongs to Fatou set ( strictly to basin of attraction of finite fixed point )
*/
// uses global variables :
// ax, ay (output = alfa(c))
double complex GiveAlfaFixedPoint(double complex c)
{
double dx, dy; //The discriminant is the d=b^2- 4ac = dx+dy*i
double r; // r(d)=sqrt(dx^2 + dy^2)
double sx, sy; // s = sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx + sy*i
double ax, ay;
// d=1-4c = dx+dy*i
dx = 1 - 4*creal(c);
dy = -4 * cimag(c);
// r(d)=sqrt(dx^2 + dy^2)
r = sqrt(dx*dx + dy*dy);
//sqrt(d) = s =sx +sy*i
sx = sqrt((r+dx)/2);
sy = sqrt((r-dx)/2);
// alfa = ax +ay*i = (1-sqrt(d))/2 = (1-sx + sy*i)/2
ax = 0.5 - sx/2.0;
ay = sy/2.0;
return ax+ay*I;
}
int setup()
{
denominator = period;
InternalAngle = 1.0/((double) denominator);
/* 2D array ranges */
if (!(iHeight % 2)) iHeight+=1; // it sholud be even number (variable % 2) or (variable & 1)
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].
iyAboveAxisLength = (iHeight -1)/2;
iyAboveMax = iyAboveAxisLength ;
iyBelowAxisLength = iyAboveAxisLength; // the same
iyAxisOfSymmetry = iyMin + iyBelowAxisLength ;
// 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))/((float)iWidth/(float)iHeight); // it should be 1.000 ...
distanceMax = PixelWidth;
// for numerical optimisation
ER2 = ER * ER;
AR = PixelWidth*15.0; // !!!! important value ( and precision and time of creation of the pgm image )
AR2= AR*AR;
PixelWidth2 = PixelWidth*PixelWidth;
/* create dynamic 1D arrays for colors ( shades of gray ) */
data = malloc( iSize * sizeof(unsigned char) );
edge = malloc( iSize * sizeof(unsigned char) );
if (edge == NULL || edge == NULL)
{
fprintf(stderr," Could not allocate memory\n");
return 1;
}
else fprintf(stderr," memory is OK \n");
return 0;
}
// from screen to world coordinate ; linear mapping
// uses global cons
double GiveZx(unsigned int ix)
{ return (ZxMin + ix*PixelWidth );}
// uses globaal cons
double GiveZy(unsigned int iy)
{ return (ZyMax - iy*PixelHeight);} // reverse y axis
// forward iteration of complex quadratic polynomial
// fc(z) = z*z +c
// initial point is a z0 = critical point
// checks if points escapes to infinity
// uses global vars : ER2, Cx, Cy, iterMax
long int GiveLastIteration(double Zx, double Zy )
{
long int iter;
// z= Zx + Zy*i
double Zx2, Zy2; /* Zx2=Zx*Zx; Zy2=Zy*Zy */
/* initial value of orbit = critical point Z= 0 */
Zx2=Zx*Zx;
Zy2=Zy*Zy;
// for iter from 0 to iterMax because of ++ after last loop
for (iter=0; iter<iterMax && ((Zx2+Zy2)<ER2); ++iter)
{
Zy=2*Zx*Zy + Cy;
Zx=Zx2-Zy2 + Cx;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
};
return iter;
}
double GiveTurn(double complex z)
{
double argument;
argument = carg(z); // argument in radians from -pi to pi
if (argument<0) argument=argument + TwoPi; // argument in radians from 0 to 2*pi
return argument/TwoPi ; // argument in turns from 0.0 to 1.0
}
/* */
int GiveNumber(double _Zx0, double _Zy0,double C_x, double C_y, int iMax, double _AR2, double _ZAx, double _ZAy )
{
int i;
double Zx, Zy; /* z = zx+zy*i */
double Zx2, Zy2; /* Zx2=Zx*Zx; Zy2=Zy*Zy */
double d, dX, dY; /* distance from z to Alpha */
double angle;
Zx=_Zx0; /* initial value of orbit */
Zy=_Zy0;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
dX=Zx-_ZAx;
dY=Zy-_ZAy;
d=dX*dX+dY*dY;
while ( d>_AR2)
{
for(i=0;i<period ;++i) // iMax = period !!!!
{
Zy=2*Zx*Zy + C_y;
Zx=Zx2-Zy2 +C_x;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
}
dX=Zx-_ZAx;
dY=Zy-_ZAy;
d=dX*dX+dY*dY;
};
// divide interior into 5 components
angle =GiveTurn(dX +dY*I) - 0.17298227404701; // this turn offset is computed with Maxima CAS filename
// star is slightly turned so this offset
// theory : Complex Dynamics by Lennart Carleson,Theodore Williams Gamelin page 40
// TurnOffset = carg(a)/q
// where a = coeff( z^(q+1))
if ((0.0<angle) && (angle<0.2)) return 0;
if ((0.2<angle) && (angle<0.4)) return 1;
if ((0.4<angle) && (angle<0.6)) return 2;
if ((0.6<angle) && (angle<0.8)) return 3;
return 4;
}
unsigned char GiveColor(unsigned int ix, unsigned int iy)
{
double Zx, Zy; // Z= Zx+ZY*i;
int iter;
unsigned char color; // gray from 0 to 255
// unsigned char ColorList[]={255,230,180}; /* shades of gray used in imaage
//color=0;
// from screen to world coordinate
Zx = GiveZx(ix);
Zy = GiveZy(iy);
iter =GiveLastIteration(Zx, Zy );
if (iter< iterMax ) // point escapes
{ color = 245; } // exterior
else // Filled Julia Set
{ int number =GiveNumber(Zx,Zy, Cx, Cy, period, AR2, creal(alfa), cimag(alfa));
color = colors[number]; }
return color;
}
/* ----------- array functions -------------- */
/* 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; }
// ix = i % iWidth;
// iy = (i- ix) / iWidth;
// i = Give_i(ix, iy);
// plots raster point (ix,iy)
int PlotPoint(unsigned int ix, unsigned int iy, unsigned char iColor)
{
unsigned i; /* index of 1D array */
i = Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */
data[i] = iColor;
return 0;
}
// fill array
// uses global var : ...
// scanning complex plane
int FillArray(unsigned char data[] )
{
unsigned int ix, iy; // pixel coordinate
// for all pixels of image
for(iy = iyMin; iy<=iyMax; ++iy)
{ printf(" %d z %d\n", iy, iyMax); //info
for(ix= ixMin; ix<=ixMax; ++ix) PlotPoint(ix, iy, GiveColor(ix, iy) ); //
}
return 0;
}
// fill array using symmetry of image
// uses global var : ...
int FillArraySymmetric(unsigned char data[] )
{
unsigned char Color; // gray from 0 to 255
printf("axis of symmetry \n");
iy = iyAxisOfSymmetry;
#pragma omp parallel for schedule(dynamic) private(ix,Color) shared(ixMin,ixMax, iyAxisOfSymmetry)
for(ix=ixMin;ix<=ixMax;++ix) {//printf(" %d from %d\n", ix, ixMax); //info
PlotPoint(ix, iy, GiveColor(ix, iy));
}
/*
The use of ‘shared(variable, variable2) specifies that these variables should be shared among all the threads.
The use of ‘private(variable, variable2)’ specifies that these variables should have a seperate instance in each thread.
*/
#pragma omp parallel for schedule(dynamic) private(iyAbove,ix,iy,Color) shared(iyAboveMin, iyAboveMax,ixMin,ixMax, iyAxisOfSymmetry)
// above and below axis
for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove)
{printf(" %d from %d\r", iyAbove, iyAboveMax); //info
for(ix=ixMin; ix<=ixMax; ++ix)
{ // above axis compute color and save it to the array
iy = iyAxisOfSymmetry + iyAbove;
Color = GiveColor(ix, iy);
PlotPoint(ix, iy, Color );
// below the axis only copy Color the same as above without computing it
PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color );
}
}
return 0;
}
int AddBoundaries(unsigned char data[])
{
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;
printf(" find boundaries in data array using Sobel filter\n");
#pragma omp parallel for schedule(dynamic) private(i,iY,iX,Gv,Gh,G) shared(iyMax,ixMax, ER2)
for(iY=1;iY<iyMax-1;++iY){
for(iX=1;iX<ixMax-1;++iX){
Gv= data[Give_i(iX-1,iY+1)] + 2*data[Give_i(iX,iY+1)] + data[Give_i(iX-1,iY+1)] - data[Give_i(iX-1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[Give_i(iX+1,iY-1)];
Gh= data[Give_i(iX+1,iY+1)] + 2*data[Give_i(iX+1,iY)] + data[Give_i(iX-1,iY-1)] - data[Give_i(iX+1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[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) {edge[i]=255;} /* background */
else {edge[i]=0;} /* boundary */
}
}
// copy boundaries from edge array to data array
for(iY=1;iY<iyMax-1;++iY){
for(iX=1;iX<ixMax-1;++iX){i= Give_i(iX,iY); if (edge[i]==0) data[i]=0;}}
return 0;
}
// Check Orientation of image : first quadrant in upper right position
// uses global var : ...
int CheckOrientation(unsigned char data[] )
{
unsigned int ix, iy; // pixel coordinate
double Zx, Zy; // Z= Zx+ZY*i;
unsigned i; /* index of 1D array */
for(iy=iyMin;iy<=iyMax;++iy)
{
Zy = GiveZy(iy);
for(ix=ixMin;ix<=ixMax;++ix)
{
// from screen to world coordinate
Zx = GiveZx(ix);
i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
if (Zx>0 && Zy>0) data[i]=255-data[i]; // check the orientation of Z-plane by marking first quadrant */
}
}
return 0;
}
// save data array to pgm file
int SaveArray2PGMFile( unsigned char data[], unsigned long int iMax)
{
FILE * fp;
const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ; it is 8 bit color file */
char name [10]; /* name of file */
sprintf(name,"new%lu", iMax); /* */
char *filename =strcat(name,".pgm");
char *comment="# ";/* comment should start with # */
/* 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(data,iSize,1,fp); /*write image data bytes to the file in one step */
printf("File %s saved. \n", filename);
fclose(fp);
return 0;
}
int info()
{
// diplay info messages
printf("Cx = %f \n", Cx);
printf("Cy = %f \n", Cy);
//
printf("alfax = %f \n", creal(alfa));
printf("alfay = %f \n", cimag(alfa));
printf("iHeght = %d \n", iHeight);
printf("distorsion of image = %f \n", ratio);
return 0;
}
/* ----------------------------------------- main -------------------------------------------------------------*/
int main()
{
// here are procedures for creating image file
setup();
c = GiveC(InternalAngle, 1.0, 1) ;
Cx=creal(c);
Cy=cimag(c);
alfa = GiveAlfaFixedPoint(c);
// compute colors of pixels = image
//FillArray( data ); // no symmetry
FillArraySymmetric(data);
AddBoundaries(data);
// CheckOrientation( data );
SaveArray2PGMFile(edge , iHeight); // save array (image) to pgm file
free(data);
free(edge);
info();
return 0;
}
New c code
/*
Adam Majewski
fraktal.republika.pl
c console progam using
* symmetry
* openMP
draw julia sets
gcc t.c -lm -Wall -fopenmp -march=native
time ./a.out
iHeight =1000
iterMax = 1000
AR = PixelWidth*15.0; // if you will increase iHeight
then increase multiplier or the time of computations will be very long
real 3m31.345s
File new2001.pgm saved.
Cx = 0.356763
Cy = 0.328582
alfax = 0.154508
alfay = 0.475528
distorsion of image = 1.000000
real 407m35.605s
*/
#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; needs also -fopenmp
/* --------------------------------- global variables and consts ------------------------------------------------------------ */
// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1
unsigned int ix, iy; // var
unsigned int ixMin = 0; // Indexes of array starts from 0 not 1
unsigned int ixMax ; //
unsigned int iWidth ; // horizontal dimension of array
unsigned int ixAxisOfSymmetry ; //
unsigned int iyMin = 0; // Indexes of array starts from 0 not 1
unsigned int iyMax ; //
unsigned int iyAxisOfSymmetry ; //
unsigned int iyAbove ; // var, measured from 1 to (iyAboveAxisLength -1)
unsigned int iyAboveMin = 1 ; //
unsigned int iyAboveMax ; //
unsigned int iyAboveAxisLength ; //
unsigned int iyBelowAxisLength ; //
unsigned int iHeight = 2000; // odd number !!!!!! = (iyMax -iyMin + 1) = iyAboveAxisLength + iyBelowAxisLength +1
// The size of array has to be a positive constant integer
unsigned int iSize ; // = iWidth*iHeight;
// memmory 1D array
unsigned char *data;
unsigned char *edge;
// unsigned int i; // var = index of 1D array
unsigned int iMin = 0; // Indexes of array starts from 0 not 1
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
/* world ( double) coordinate = dynamic plane */
const double ZxMin=-1.5;
const double ZxMax=1.5;
const double ZyMin=-1.5;
const double ZyMax=1.5;
double PixelWidth; // =(ZxMax-ZxMin)/iXmax;
double PixelWidth2; // = PixelWidth*PixelWidth;
double PixelHeight; // =(ZyMax-ZyMin)/iYmax;
double distanceMax;
double ratio ;
double TwoPi=2*M_PI;
// complex numbers of parametr plane
double Cx; // c =Cx +Cy * i
double Cy;
double complex c; //
double complex alfa; // alfa fixed point alfa=f(alfa)
unsigned long int iterMax = 1000; //iHeight*100;
double ER = 2.0; // Escape Radius for bailout test
double ER2;
double AR ; // radius around attractor
double AR2; // =AR*AR;
unsigned int period = 5; // period of secondary component joined by root point
unsigned int denominator;
double InternalAngle;
unsigned char colors[]={255,230,180, 160,140}; // number of colors >= period
/* ------------------------------------------ functions -------------------------------------------------------------*/
/* find c in component of Mandelbrot set
uses code by Wolf Jung from program Mandel
see function mndlbrot::bifurcate from mandelbrot.cpp
http://www.mndynamics.com/indexp.html
*/
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int period)
{
//0 <= InternalRay<= 1
//0 <= InternalAngleInTurns <=1
double t = InternalAngleInTurns *2*M_PI; // from turns to radians
double R2 = InternalRadius * InternalRadius;
//double Cx, Cy; /* C = Cx+Cy*i */
switch ( period ) // of component
{
case 1: // main cardioid
Cx = (cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4;
Cy = (sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4;
break;
case 2: // only one component
Cx = InternalRadius * 0.25*cos(t) - 1.0;
Cy = InternalRadius * 0.25*sin(t);
break;
// for each period there are 2^(period-1) roots.
default: // higher periods : to do
Cx = 0.0;
Cy = 0.0;
break; }
return Cx + Cy*I;
}
/*
http://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings
z^2 + c = z
z^2 - z + c = 0
ax^2 +bx + c =0 // ge3neral for of quadratic equation
so :
a=1
b =-1
c = c
so :
The discriminant is the d=b^2- 4ac
d=1-4c = dx+dy*i
r(d)=sqrt(dx^2 + dy^2)
sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx +- sy*i
x1=(1+sqrt(d))/2 = beta = (1+sx+sy*i)/2
x2=(1-sqrt(d))/2 = alfa = (1-sx -sy*i)/2
alfa : attracting when c is in main cardioid of Mandelbrot set, then it is in interior of Filled-in Julia set,
it means belongs to Fatou set ( strictly to basin of attraction of finite fixed point )
*/
// uses global variables :
// ax, ay (output = alfa(c))
double complex GiveAlfaFixedPoint(double complex c)
{
double dx, dy; //The discriminant is the d=b^2- 4ac = dx+dy*i
double r; // r(d)=sqrt(dx^2 + dy^2)
double sx, sy; // s = sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx + sy*i
double ax, ay;
// d=1-4c = dx+dy*i
dx = 1 - 4*creal(c);
dy = -4 * cimag(c);
// r(d)=sqrt(dx^2 + dy^2)
r = sqrt(dx*dx + dy*dy);
//sqrt(d) = s =sx +sy*i
sx = sqrt((r+dx)/2);
sy = sqrt((r-dx)/2);
// alfa = ax +ay*i = (1-sqrt(d))/2 = (1-sx + sy*i)/2
ax = 0.5 - sx/2.0;
ay = sy/2.0;
return ax+ay*I;
}
int setup()
{
denominator = period;
InternalAngle = 1.0/((double) denominator);
/* 2D array ranges */
if (!(iHeight % 2)) iHeight+=1; // it sholud be even number (variable % 2) or (variable & 1)
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].
iyAboveAxisLength = (iHeight -1)/2;
iyAboveMax = iyAboveAxisLength ;
iyBelowAxisLength = iyAboveAxisLength; // the same
iyAxisOfSymmetry = iyMin + iyBelowAxisLength ;
// 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))/((float)iWidth/(float)iHeight); // it should be 1.000 ...
distanceMax = PixelWidth;
// for numerical optimisation
ER2 = ER * ER;
AR = PixelWidth*50.0; // !!!! important value ( and precision and time of creation of the pgm image )
AR2= AR*AR;
PixelWidth2 = PixelWidth*PixelWidth;
/* create dynamic 1D arrays for colors ( shades of gray ) */
data = malloc( iSize * sizeof(unsigned char) );
edge = malloc( iSize * sizeof(unsigned char) );
if (edge == NULL || edge == NULL)
{
fprintf(stderr," Could not allocate memory\n");
return 1;
}
else fprintf(stderr," memory is OK \n");
return 0;
}
// from screen to world coordinate ; linear mapping
// uses global cons
double GiveZx(unsigned int ix)
{ return (ZxMin + ix*PixelWidth );}
// uses globaal cons
double GiveZy(unsigned int iy)
{ return (ZyMax - iy*PixelHeight);} // reverse y axis
// forward iteration of complex quadratic polynomial
// fc(z) = z*z +c
// initial point is a z0 = critical point
// checks if points escapes to infinity
// uses global vars : ER2, Cx, Cy, iterMax
long int GiveLastIteration(double Zx, double Zy )
{
long int iter;
// z= Zx + Zy*i
double Zx2, Zy2; /* Zx2=Zx*Zx; Zy2=Zy*Zy */
/* initial value of orbit = critical point Z= 0 */
Zx2=Zx*Zx;
Zy2=Zy*Zy;
// for iter from 0 to iterMax because of ++ after last loop
for (iter=0; iter<iterMax && ((Zx2+Zy2)<ER2); ++iter)
{
Zy=2*Zx*Zy + Cy;
Zx=Zx2-Zy2 + Cx;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
};
return iter;
}
double GiveTurn(double complex z)
{
double argument;
argument = carg(z); // argument in radians from -pi to pi
if (argument<0) argument=argument + TwoPi; // argument in radians from 0 to 2*pi
return argument/TwoPi ; // argument in turns from 0.0 to 1.0
}
/* */
int GiveNumber(double _Zx0, double _Zy0,double C_x, double C_y, int iMax, double _AR2, double _ZAx, double _ZAy )
{
int i;
double Zx, Zy; /* z = zx+zy*i */
double Zx2, Zy2; /* Zx2=Zx*Zx; Zy2=Zy*Zy */
double d, dX, dY; /* distance from z to Alpha */
double angle;
Zx=_Zx0; /* initial value of orbit */
Zy=_Zy0;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
dX=Zx-_ZAx;
dY=Zy-_ZAy;
d=dX*dX+dY*dY;
while ( d>_AR2)
{
for(i=0;i<period ;++i) // iMax = period !!!!
{
Zy=2*Zx*Zy + C_y;
Zx=Zx2-Zy2 +C_x;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
}
dX=Zx-_ZAx;
dY=Zy-_ZAy;
d=dX*dX+dY*dY;
};
// divide interior into 5 components
angle =GiveTurn(dX +dY*I); // - 0.17298227404701; // this turn offset is computed with Maxima CAS filename
// star is slightly turned so this offset
// theory : Complex Dynamics by Lennart Carleson,Theodore Williams Gamelin page 40
// TurnOffset = carg(a)/q
// where a = coeff( z^(q+1))
if (angle < 0.16 || angle > 0.97) return 0; // 0.0 < angle < 0.19
if (angle < 0.36 ) return 1; // 0.19 < angle
if (angle < 0.57 ) return 2; // ((0.4<angle) &&
if (angle < 0.775 ) return 3; // ((0.6<angle) &&
return 4;
}
unsigned char GiveColor(unsigned int ix, unsigned int iy)
{
double Zx, Zy; // Z= Zx+ZY*i;
int iter;
unsigned char color; // gray from 0 to 255
// unsigned char ColorList[]={255,230,180}; /* shades of gray used in imaage
//color=0;
// from screen to world coordinate
Zx = GiveZx(ix);
Zy = GiveZy(iy);
iter =GiveLastIteration(Zx, Zy );
if (iter< iterMax ) // point escapes
{ color = 245; } // exterior
else // Filled Julia Set
{ int number =GiveNumber(Zx,Zy, Cx, Cy, period, AR2, creal(alfa), cimag(alfa));
color = colors[number]; }
return color;
}
/* ----------- array functions -------------- */
/* 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; }
// ix = i % iWidth;
// iy = (i- ix) / iWidth;
// i = Give_i(ix, iy);
// plots raster point (ix,iy)
int PlotPoint(unsigned int ix, unsigned int iy, unsigned char iColor)
{
unsigned i; /* index of 1D array */
i = Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */
data[i] = iColor;
return 0;
}
// fill array
// uses global var : ...
// scanning complex plane
int FillArray(unsigned char data[] )
{
unsigned int ix, iy; // pixel coordinate
// for all pixels of image
for(iy = iyMin; iy<=iyMax; ++iy)
{ printf(" %d z %d\n", iy, iyMax); //info
for(ix= ixMin; ix<=ixMax; ++ix) PlotPoint(ix, iy, GiveColor(ix, iy) ); //
}
return 0;
}
// fill array using symmetry of image
// uses global var : ...
int FillArraySymmetric(unsigned char data[] )
{
unsigned char Color; // gray from 0 to 255
printf("axis of symmetry \n");
iy = iyAxisOfSymmetry;
#pragma omp parallel for schedule(dynamic) private(ix,Color) shared(ixMin,ixMax, iyAxisOfSymmetry)
for(ix=ixMin;ix<=ixMax;++ix) {//printf(" %d from %d\n", ix, ixMax); //info
PlotPoint(ix, iy, GiveColor(ix, iy));
}
/*
The use of ‘shared(variable, variable2) specifies that these variables should be shared among all the threads.
The use of ‘private(variable, variable2)’ specifies that these variables should have a seperate instance in each thread.
*/
#pragma omp parallel for schedule(dynamic) private(iyAbove,ix,iy,Color) shared(iyAboveMin, iyAboveMax,ixMin,ixMax, iyAxisOfSymmetry)
// above and below axis
for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove)
{printf(" %d from %d\r", iyAbove, iyAboveMax); //info
for(ix=ixMin; ix<=ixMax; ++ix)
{ // above axis compute color and save it to the array
iy = iyAxisOfSymmetry + iyAbove;
Color = GiveColor(ix, iy);
PlotPoint(ix, iy, Color );
// below the axis only copy Color the same as above without computing it
PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color );
}
}
return 0;
}
int AddBoundaries(unsigned char data[])
{
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;
printf(" find boundaries in data array using Sobel filter\n");
#pragma omp parallel for schedule(dynamic) private(i,iY,iX,Gv,Gh,G) shared(iyMax,ixMax, ER2)
for(iY=1;iY<iyMax-1;++iY){
for(iX=1;iX<ixMax-1;++iX){
Gv= data[Give_i(iX-1,iY+1)] + 2*data[Give_i(iX,iY+1)] + data[Give_i(iX-1,iY+1)] - data[Give_i(iX-1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[Give_i(iX+1,iY-1)];
Gh= data[Give_i(iX+1,iY+1)] + 2*data[Give_i(iX+1,iY)] + data[Give_i(iX-1,iY-1)] - data[Give_i(iX+1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[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) {edge[i]=255;} /* background */
else {edge[i]=0;} /* boundary */
}
}
// copy boundaries from edge array to data array
for(iY=1;iY<iyMax-1;++iY){
for(iX=1;iX<ixMax-1;++iX){i= Give_i(iX,iY); if (edge[i]==0) data[i]=0;}}
return 0;
}
// Check Orientation of image : first quadrant in upper right position
// uses global var : ...
int CheckOrientation(unsigned char data[] )
{
unsigned int ix, iy; // pixel coordinate
double Zx, Zy; // Z= Zx+ZY*i;
unsigned i; /* index of 1D array */
for(iy=iyMin;iy<=iyMax;++iy)
{
Zy = GiveZy(iy);
for(ix=ixMin;ix<=ixMax;++ix)
{
// from screen to world coordinate
Zx = GiveZx(ix);
i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
if (Zx>0 && Zy>0) data[i]=255-data[i]; // check the orientation of Z-plane by marking first quadrant */
}
}
return 0;
}
// save data array to pgm file
int SaveArray2PGMFile( unsigned char data[], unsigned long int iMax)
{
FILE * fp;
const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ; it is 8 bit color file */
char name [10]; /* name of file */
sprintf(name,"new%lu", iMax); /* */
char *filename =strcat(name,".pgm");
char *comment="# ";/* comment should start with # */
/* 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(data,iSize,1,fp); /*write image data bytes to the file in one step */
printf("File %s saved. \n", filename);
fclose(fp);
return 0;
}
int info()
{
// diplay info messages
printf("Cx = %f \n", Cx);
printf("Cy = %f \n", Cy);
//
printf("alfax = %f \n", creal(alfa));
printf("alfay = %f \n", cimag(alfa));
printf("iHeght = %d \n", iHeight);
printf("distorsion of image = %f \n", ratio);
return 0;
}
/* ----------------------------------------- main -------------------------------------------------------------*/
int main()
{
// here are procedures for creating image file
setup();
c = GiveC(InternalAngle, 1.0, 1) ;
Cx=creal(c);
Cy=cimag(c);
alfa = GiveAlfaFixedPoint(c);
// compute colors of pixels = image
//FillArray( data ); // no symmetry
FillArraySymmetric(data);
AddBoundaries(data);
// CheckOrientation( data );
SaveArray2PGMFile(edge , iHeight); // save array (image) to pgm file
free(data);
free(edge);
info();
return 0;
}
result of new code :
memory is OK axis of symmetry find boundaries in data array using Sobel filter File new2001.pgm saved. Cx = 0.356763 Cy = 0.328582 alfax = 0.154508 alfay = 0.475528 iHeght = 2001 distorsion of image = 1.000000 real 0m7.806s user 1m1.933s sys 0m0.042s
Rererences
- ↑ Complex Dynamics by Lennart Carleson,Theodore Williams Gamelin page 40
File history
Click on a date/time to view the file as it appeared at that time.
| Date/Time | Thumbnail | Dimensions | User | Comment | |
|---|---|---|---|---|---|
| current | 17:28, 27 December 2012 | | 2,001 × 2,001 (33 KB) | Soul windsurfer | {{Information |Description ={{en|1=Parabolic Julia set for internal angle 1 over 5 - new method}} |Source ={{own}} |Author =Adam majewski |Date =2012-12-26 |Permission = |other_versions = }} |
File usage
The following 2 pages use this file:
