File:Quadratic Julia set with Internal binary decomposition for internal ray 1 over 3.ogv
Jump to navigation
Jump to search
Size of this JPG preview of this OGG file: 600 × 600 pixels. Other resolutions: 240 × 240 pixels | 480 × 480 pixels | 1,001 × 1,001 pixels.
Original file (Ogg Theora video file, length 12 s, 1,001 × 1,001 pixels, 3.8 Mbps, file size: 5.34 MB)
Captions
Captions
Add a one-line explanation of what this file represents
Summary[edit]
| Description |
English: Quadratic Julia set with Internal binary decomposition for internal ray 1 over 3 |
| Date | |
| Source | Own work |
| Author | Adam majewski |
Summary[edit]
This video is inspired by images by T Kawahira [1]
Meaning of colors of interior of filled Julia set:
- " distinguish the regions mapped to distinct copies of C in the linearized models (§4)"[2]
- show which points form components of Julia set as c moves along internal ray 1/3. See video
Number in the left upper corner is a internal radius.
Algorithm[edit]
- scan all points of dynamical plane ( see function FillArraySymmetric ). Use symmetry of image and all cores of CPU ( OpenMP ) to speed up computations.
- For each point of dynamical plane ( = memory array data ) compute it's colour and save it to memory array ( see function GiveColor).
- add boundaries to image ( see function AddBoundaries )
- Save array to graphic file ( see function SaveArray2PGMFile )
Colour[edit]
- First divide all points into 2 sets ( using bailout test = function GiveLastIteration ) :
- basin of attraction to infinity ( exterior of filled Julia set ) = all point that escapes
- filled Julia set = all points that do not escapes = bounded by circle with radius ER : abs(z)<ER
- All points that do not escapes treat as points of interior of filled Julia set which tend to alfa fixed point. Divide them into 3 subsets ( see function GiveNumber ):
- 0 < turn of zn < 1/3
- 1/3 <= turn of zn < 2/3
- 2/3 < turn of zn <=1
// divide interior into 3 components
angle =GiveTurn(dX +dY*I);
if ((0.0<angle) && (angle<0.333333333)) return 0;
if ((0.333333333<angle) && (angle<0.666666667)) return 1;
return 2;
Turn of zn =
The trick is to check turn of zn :
- when zn is inside target set ( trap = set around alfa fixed point ZA with radius AR )
- after n*period iterations where n is positive integer and period denotes period of next hyperbolic component of Mandelbrot set to which internal ray goes.
In this case angle of internal ray is 1/3 and it goes to period 3 hyperbolic component. Not that c i inside main cardioid of Mandelbrot set ( period = 1 )
while ( d>_AR2)
{
for(i=0;i<3;++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;
};
C source code[edit]
This C code creates iMax pgm files :
/*
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
components of interior are invariant under f^period
*/
#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 = 1000; // 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 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 = 3; // period of secondary component joined by root point
unsigned int denominator;
double InternalAngle;
unsigned char colors[]={255,230,180};
/* ------------------------------------------ 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;
PixelWidth2 = PixelWidth*PixelWidth;
/* create dynamic 1D arrays for colors ( shades of gray ) */
data = malloc( iSize * sizeof(unsigned char) );
if (data == NULL )
{
fprintf(stderr," Could not allocate memory\n");
getchar();
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<3;++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 3 components
angle =GiveTurn(dX +dY*I);
if ((0.0<angle) && (angle<0.333333333)) return 0;
if ((0.333333333<angle) && (angle<0.666666667)) return 1;
return 2;
}
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;
// memmory 1D array
unsigned char *edge;
/* create dynamic 1D arrays for colors ( shades of gray ) */
edge = malloc( iSize * sizeof(unsigned char) );
if (edge == NULL )
{
fprintf(stderr," Could not allocate memory\n");
return 1;
}
// 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 */
}
}
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;}}
free(edge);
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[],double Radius)
{
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,"%f", Radius); /* */
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("distorsion of image = %f \n", ratio);
return 0;
}
double GiveAR2(double R, int i, int iMax, double PixelWidth )
{
if (R>=0.09)
AR = PixelWidth*15;// AR is const
else AR = PixelWidth*(i/(0.06*iMax)+1);// AR is related with internal radius
return (AR*AR);
}
/* ----------------------------------------- main -------------------------------------------------------------*/
int main()
{
// here are procedures for creating image file
int i;
int iMax=100;
double stepR =1.0/iMax;
double R= 0.0;
setup();
for(i=0;i<=iMax;++i)
{
//printf("R= %f \n", R);
R=i*stepR;
c = GiveC(InternalAngle, R, 1) ;
Cx=creal(c);
Cy=cimag(c);
alfa = GiveAlfaFixedPoint(c);
// AR = PixelWidth*15;// AR is const
AR2= GiveAR2(R, i, iMax,PixelWidth);
// compute colors of pixels = image
//FillArray( data ); // no symmetry
FillArraySymmetric(data);
AddBoundaries(data);
// CheckOrientation( data );
SaveArray2PGMFile( data,R); // save array (image) to pgm file
}
free(data);
//info();
return 0;
}
Bash source code[edit]
#!/bin/bash
# script file for BASH
# which bash
# save this file as g
# chmod +x g
# ./g
i=0
# for all pgm files in this directory
for file in *.pgm ; do
# b is name of file without extension
b=$(basename $file .pgm)
# change file name to integers and count files
((i= i+1))
# convert from pgm to gif and add text ( level ) using ImageMagic
convert $file -pointsize 50 -annotate +10+100 $b ${i}.gif
echo $file
done
echo convert all gif files to one video file
# ffmpeg2theora %d.gif --framerate 5 --videoquality 9 -f webm --artist "Adam Majewski" -o o${i}.webm
ffmpeg2theora %d.gif --framerate 10 --videoquality 9 -f ogv --artist "Adam Majewski" -o o${i}.ogv
echo o${i} OK
# end
References[edit]
- ↑ images by T Kawahira[dead link]
- ↑ Semiconjugacies in Complex Dynamics with Parabolics by Tomoki Kawahira. ( page 49 ) archive copy at the Wayback Machine
Licensing[edit]
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.
File history
Click on a date/time to view the file as it appeared at that time.
| Date/Time | Thumbnail | Dimensions | User | Comment | |
|---|---|---|---|---|---|
| current | 13:40, 24 December 2012 | 12 s, 1,001 × 1,001 (5.34 MB) | Soul windsurfer (talk | contribs) | better quality | |
| 18:59, 23 December 2012 | 12 s, 1,001 × 1,001 (3.74 MB) | Soul windsurfer (talk | contribs) | {{Information |Description ={{en|1=Quadratic Julia set with Internal binary decomposition for internal ray 1 over 3}} |Source ={{own}} |Author =Adam majewski |Date =2012-12-23 |Permission = |other... |
You cannot overwrite this file.
File usage on Commons
The following page uses this file:
Transcode status
Update transcode statusFile usage on other wikis
The following other wikis use this file:
- Usage on en.wikibooks.org