Plik:Escher like tiling Fat Douady rabbit Julia set.png
Przejdź do nawigacji
Przejdź do wyszukiwania
Rozmiar podglądu – 600 × 600 pikseli. Inne rozdzielczości: 240 × 240 pikseli | 480 × 480 pikseli | 768 × 768 pikseli | 1024 × 1024 pikseli | 2000 × 2000 pikseli.
Rozmiar pierwotny (2000 × 2000 pikseli, rozmiar pliku: 646 KB, typ MIME: image/png)
Opis
| Opis |
English: Escher like tiling Fat Douady rabbit Julia set. c = -0.125+0.649519052838329*i. Compare Fig 4.12 page 185 from book : The Science of Fractal Images by MF Barnsley et al.[1] |
| Data | |
| Źródło | own work with help of Michael Sargent |
| Autor | Adam majewski |
| Inne wersje |
|
Licencja
Ja, właściciel praw autorskich do tego dzieła, udostępniam je na poniższej licencji
Ten plik udostępniony jest na licencji Creative Commons Uznanie autorstwa – Na tych samych warunkach 4.0 Międzynarodowe.
- Wolno:
- dzielić się – kopiować, rozpowszechniać, odtwarzać i wykonywać utwór
- modyfikować – tworzyć utwory zależne
- Na następujących warunkach:
- uznanie autorstwa – musisz określić autorstwo utworu, podać link do licencji, a także wskazać czy utwór został zmieniony. Możesz to zrobić w każdy rozsądny sposób, o ile nie będzie to sugerować, że licencjodawca popiera Ciebie lub Twoje użycie utworu.
- na tych samych warunkach – Jeśli zmienia się lub przekształca niniejszy utwór, lub tworzy inny na jego podstawie, można rozpowszechniać powstały w ten sposób nowy utwór tylko na podstawie tej samej lub podobnej licencji.
c source code
/*
Adam Majewski
adammaj1 aaattt o2 dot pl // o like oxygen not 0 like zero
==============================================
Structure of a program or how to analyze the program
============== Image X ========================
DrawImageOfX -> DrawPointOfX -> ComputeColorOfX
first 2 functions are identical for every X
check only last function = ComputeColorOfX
which computes color of one pixel !
==========================================
---------------------------------
indent d.c
default is gnu style
-------------------
c console progam
export OMP_DISPLAY_ENV="TRUE"
gcc d.c -lm -Wall -march=native -fopenmp
time ./a.out > b.txt
gcc d.c -lm -Wall -march=native -fopenmp
time ./a.out
time ./a.out >a.txt
----------------------
real 0m19,809s
user 2m26,763s
sys 0m0,161s
{c=0+0*i; circle }
C_re=-0.0; {wspolczynnik funkcji F(z)=(z*z)+c }
C_im=0.0; { c jest liczba zespolona : c=c_re+c_im*i }
{target set C_int }
C_re_int=-1.1; {wspolczynnik funkcji F(z)=(z*z)+c_int }
C_im_int=-0.001; { c jest liczba zespolona : c_int=c_re_int+c_im_int*i }
skala=250.10; {standard 15 ; sproboj: -2.0; 1.0 ; 1000.0; }
xMin=-1.0;
xMax=1.0;
dx=xMax-xMin;
yMin=-1.0;
yMax=1.0;
*/
#include <stdio.h>
#include <stdlib.h> // malloc
#include <string.h> // strcat
#include <math.h> // M_PI; needs -lm also
#include <complex.h>
#include <omp.h> // OpenMP
/* --------------------------------- global variables and consts ------------------------------------------------------------ */
// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1
//unsigned int ix, iy; // var
static unsigned int ixMin = 0; // Indexes of array starts from 0 not 1
static unsigned int ixMax; //
static unsigned int iWidth; // horizontal dimension of array
static unsigned int iyMin = 0; // Indexes of array starts from 0 not 1
static unsigned int iyMax; //
static unsigned int iHeight = 10000; //
// The size of array has to be a positive constant integer
static unsigned int iSize; // = iWidth*iHeight;
// memmory 1D array
unsigned char *data;
unsigned char *edge;
//unsigned char *edge2;
// unsigned int i; // var = index of 1D array
//static unsigned int iMin = 0; // Indexes of array starts from 0 not 1
static unsigned int iMax; // = i2Dsize-1 =
// The size of array has to be a positive constant integer
// unsigned int i1Dsize ; // = i2Dsize = (iMax -iMin + 1) = ; 1D array with the same size as 2D array
static const double ZxMin = -1.0; //-0.05;
static const double ZxMax = 1.0; //0.75;
static const double ZyMin = -1.0; //-0.1;
static const double ZyMax = 1.0; //0.7;
static double PixelWidth; // =(ZxMax-ZxMin)/ixMax;
static double PixelHeight; // =(ZyMax-ZyMin)/iyMax;
static double ratio;
// complex numbers of parametr plane
double complex c = 0.0; // parameter of function fc(z)=z^2 + c
//target set C_int }
double complex tc =-0.125+0.649519052838329*I; // {wspolczynnik funkcji F(z)=(z*z)+c_int }
double scale = 10.0; // {standard 15 ; sproboj: -2.0; 1.0 ; 1000.0; }
int Period = 2;
static int iterMax = 255; //iHeight*100;
static double ER = 2.0; // EscapeRadius for bailout test
//double EscapeRadius=1000000; // = ER big !!!!
// SAC/J
//double lnER; // ln(ER)
//int i_skip = 2; // exclude (i_skip+1) elements from average
//unsigned char s = 7; // stripe density
//double BoundaryWidth = 3.0; // % of image width
//double distanceMax; //distanceMax = BoundaryWidth*PixelWidth;
/* colors = shades of gray from 0 to 255 */
unsigned char iColorOfExterior = 250;
unsigned char iColorOfInterior = 200;
unsigned char iColorOfInterior1 = 210;
unsigned char iColorOfInterior2 = 180;
unsigned char iColorOfBoundary = 0;
unsigned char iColorOfUnknown = 30;
/* ------------------------------------------ functions -------------------------------------------------------------*/
//------------------complex numbers -----------------------------------------------------
// from screen to world coordinate ; linear mapping
// uses global cons
double GiveZx ( int ix)
{
return (ZxMin + ix * PixelWidth);
}
// uses globaal cons
double GiveZy (int iy) {
return (ZyMax - iy * PixelHeight);
} // reverse y axis
complex double GiveZ( int ix, int iy){
double Zx = GiveZx(ix);
double Zy = GiveZy(iy);
return Zx + Zy*I;
}
// ****************** DYNAMICS = trap tests ( target sets) ****************************
// bailout test
// z escapes when
// abs(z)> ER or cabs2(z)> ER2
// https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Julia_set#Boolean_Escape_time
int Escapes(complex double z){
// here target set (trap) is the exterior circle with radsius = ER ( EscapeRadius)
// with ceter = origin z= 0
// on the Riemann sphere it is a circle with point at infinity as a center
if (cabs(z)>ER) return 1;
return 0;
}
/* ----------- array functions = drawing -------------- */
/* gives position of 2D point (ix,iy) in 1D array ; uses also global variable iWidth */
unsigned int Give_i (unsigned int ix, unsigned int iy)
{
return ix + iy * iWidth;
}
// ***********************************************************************************************
// ********************** edge detection usung Sobel filter ***************************************
// ***************************************************************************************************
// from Source to Destination
int ComputeBoundaries(unsigned char S[], unsigned char D[])
{
unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
unsigned int i; /* index of 1D array */
/* sobel filter */
unsigned char G, Gh, Gv;
// boundaries are in D array ( global var )
// clear D array
memset(D, iColorOfExterior, iSize*sizeof(*D)); // for heap-allocated arrays, where N is the number of elements = FillArrayWithColor(D , iColorOfExterior);
// printf(" find boundaries in S array using Sobel filter\n");
#pragma omp parallel for schedule(dynamic) private(i,iY,iX,Gv,Gh,G) shared(iyMax,ixMax)
for(iY=1;iY<iyMax-1;++iY){
for(iX=1;iX<ixMax-1;++iX){
Gv= S[Give_i(iX-1,iY+1)] + 2*S[Give_i(iX,iY+1)] + S[Give_i(iX-1,iY+1)] - S[Give_i(iX-1,iY-1)] - 2*S[Give_i(iX-1,iY)] - S[Give_i(iX+1,iY-1)];
Gh= S[Give_i(iX+1,iY+1)] + 2*S[Give_i(iX+1,iY)] + S[Give_i(iX-1,iY-1)] - S[Give_i(iX+1,iY-1)] - 2*S[Give_i(iX-1,iY)] - S[Give_i(iX-1,iY-1)];
G = sqrt(Gh*Gh + Gv*Gv);
i= Give_i(iX,iY); /* compute index of 1D array from indices of 2D array */
if (G==0) {D[i]=255;} /* background */
else {D[i]=0;} /* boundary */
}
}
return 0;
}
// copy from Source to Destination
int CopyBoundaries(unsigned char S[], unsigned char D[])
{
unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
unsigned int i; /* index of 1D array */
//printf("copy boundaries from S array to D array \n");
for(iY=1;iY<iyMax-1;++iY)
for(iX=1;iX<ixMax-1;++iX)
{i= Give_i(iX,iY); if (S[i]==0) D[i]=0;}
return 0;
}
// ***************************************************************************************************************************
// ************************** LSM/J*****************************************
// ****************************************************************************************************************************
/*
Begin { petla wewnetrzna}
tn:=0;
tx:=x*skala;
ty:=y*skala;
repeat
inc(tn);
xTemp:=(tx+ty)*(tx-ty)+C_re_int;
ty:=2*tx*ty+C_im_int;
tx:=xTemp;
until (tn>=nMax) or (ty*ty+tx*tx>BailOut);
if tn=nMax then break;
end;
*/
int EscapeTimeEscher(complex double z){
int nMax = iterMax;
int n;
//double cabsz;
z = z*scale; //
for (n=0; n <= nMax && cabs(z)<ER; n++){ //forward iteration
//if (cabsz > ER) break; // escaping
//if (cabsz< PixelWidth) break; // fails into finite attractor = interior
z = z*z + tc; /* forward iteration : reversed basilica */
}
return n;
}
unsigned char ComputeColorOfLSM(complex double z){
int nMax = iterMax;
//double cabsz;
unsigned char iColor;
complex double tz;
int n = 0;
int tn = 0;
double size = 0.0;
double tsize = 0.0;
size = cabs(z);
while (n < nMax && (size < ER)){ //forward iteration
// inner loop
tn = 0;
tsize = 0.0;
tz = z*scale;
while (tn < nMax && (tsize < ER)){ //forward iteration
tz = tz*tz + tc; /* forward iteration complex quadratic polynomial */
tsize = cabs(tz);
tn++;
}
if (tn == nMax) break;
//
z = z*z; /* forward iteration complex quadratic polynomial */
size = cabs(z);
n++;
}
if (tn == nMax) // inside circle
{
if (n==0)
{iColor = 255;} // target set
else iColor = (n*30) % 255;
}
else iColor = iColorOfExterior; // outside the circle
return iColor;
}
// plots raster point (ix,iy)
int DrawPointOfLSM (unsigned char A[], int ix, int iy)
{
int i; /* index of 1D array */
unsigned char iColor;
complex double z;
i = Give_i (ix, iy); /* compute index of 1D array from indices of 2D array */
z = GiveZ(ix,iy);
iColor = ComputeColorOfLSM(z);
A[i] = iColor ; // interior
return 0;
}
// fill array
// uses global var : ...
// scanning complex plane
int DrawImagerOfLSM (unsigned char A[])
{
unsigned int ix, iy; // pixel coordinate
//printf("compute image \n");
// for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
for (iy = iyMin; iy <= iyMax; ++iy){
printf (" %d from %d \r", iy, iyMax); //info
for (ix = ixMin; ix <= ixMax; ++ix)
DrawPointOfLSM(A, ix, iy); //
}
return 0;
}
// ***************************************************************************************************************************
// ************************** binary decomposition BD/J*****************************************
// ****************************************************************************************************************************
unsigned char ComputeColorOfBD(complex double z){
int nMax = 255;
double cabsz;
unsigned char iColor;
int n;
for (n=0; n < nMax; n++){ //forward iteration
cabsz = cabs(z);
if (cabsz > ER) break; // esacping
if (cabsz< PixelWidth) break; // fails into finite attractor = interior
z = z*z +c ; /* forward iteration : complex quadratic polynomial */
}
if (creal(z)>0.0)
iColor = 255;
else iColor = 0;
return iColor;
}
// plots raster point (ix,iy)
int DrawPointOfBD (unsigned char A[], int ix, int iy)
{
int i; /* index of 1D array */
unsigned char iColor;
complex double z;
i = Give_i (ix, iy); /* compute index of 1D array from indices of 2D array */
z = GiveZ(ix,iy);
iColor = ComputeColorOfBD(z);
A[i] = iColor ; // interior
return 0;
}
// fill array
// uses global var : ...
// scanning complex plane
int DrawImagerOfBD (unsigned char A[])
{
unsigned int ix, iy; // pixel coordinate
//printf("compute image \n");
// for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
for (iy = iyMin; iy <= iyMax; ++iy){
printf (" %d from %d \r", iy, iyMax); //info
for (ix = ixMin; ix <= ixMax; ++ix)
DrawPointOfBD(A, ix, iy); //
}
return 0;
}
// *******************************************************************************************
// ********************************** save A array to pgm file ****************************
// *********************************************************************************************
int SaveArray2PGMFile( unsigned char A[], double k, char* comment )
{
FILE * fp;
const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ; it is 8 bit color file */
char name [100]; /* name of file */
snprintf(name, sizeof name, "%.1f", k); /* */
char *filename =strncat(name,".pgm", 4);
// save image to the pgm file
fp= fopen(filename,"wb"); // create new file,give it a name and open it in binary mode
fprintf(fp,"P5\n # %s\n %u %u\n %u\n", comment, iWidth, iHeight, MaxColorComponentValue); // write header to the file
fwrite(A,iSize,1,fp); // write array with image data bytes to the file in one step
fclose(fp);
// info
printf("File %s saved ", filename);
if (comment == NULL || strlen(comment) ==0)
printf("\n");
else printf (". Comment = %s \n", comment);
return 0;
}
int PrintInfoAboutProgam()
{
// display info messages
printf ("Numerical approximation of Escher like tilings of Julia sets by Michael Sargent for fc(z)= z^2 + c \n");
//printf ("iPeriodParent = %d \n", iPeriodParent);
//printf ("iPeriodOfChild = %d \n", iPeriodChild);
printf ("parameter c = ( %.16f ; %.16f ) \n", creal(tc), cimag(tc));
printf ("Image Width = %f in world coordinate\n", ZxMax - ZxMin);
printf ("PixelWidth = %f \n", PixelWidth);
// image corners in world coordinate
// center and radius
// center and zoom
// GradientRepetition
printf ("Maximal number of iterations = iterMax = %d \n", iterMax);
printf ("ratio of image = %f ; it should be 1.000 ...\n", ratio);
//
printf("gcc version: %d.%d.%d\n",__GNUC__,__GNUC_MINOR__,__GNUC_PATCHLEVEL__); // https://stackoverflow.com/questions/20389193/how-do-i-check-my-gcc-c-compiler-version-for-my-eclipse
// OpenMP version is diplayed in the console
return 0;
}
// *****************************************************************************
//;;;;;;;;;;;;;;;;;;;;;; setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
// **************************************************************************************
int setup ()
{
printf ("setup start\n");
c = 0.0; //
/* 2D array ranges */
iWidth = iHeight;
iSize = iWidth * iHeight; // size = number of points in array
// iy
iyMax = iHeight - 1; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].
//ix
ixMax = iWidth - 1;
/* 1D array ranges */
// i1Dsize = i2Dsize; // 1D array with the same size as 2D array
iMax = iSize - 1; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].
/* Pixel sizes */
PixelWidth = (ZxMax - ZxMin) / ixMax; // ixMax = (iWidth-1) step between pixels in world coordinate
PixelHeight = (ZyMax - ZyMin) / iyMax;
ratio = ((ZxMax - ZxMin) / (ZyMax - ZyMin)) / ((double) iWidth / (double) iHeight); // it should be 1.000 ...
//ER2 = ER * ER; // for numerical optimisation in iteration
//lnER = log(EscapeRadius); // ln(ER)
/* create dynamic 1D arrays for colors ( shades of gray ) */
data = malloc (iSize * sizeof (unsigned char));
edge = malloc (iSize * sizeof (unsigned char));
//edge2 = malloc (iSize * sizeof (unsigned char));
if (data == NULL || edge == NULL ){
fprintf (stderr, " Could not allocate memory");
return 1;
}
//BoundaryWidth = 20.0*iWidth/2000.0 ; // measured in pixels ( when iWidth = 2000)
//distanceMax = BoundaryWidth*PixelWidth;
printf (" end of setup \n");
return 0;
} // ;;;;;;;;;;;;;;;;;;;;;;;;; end of the setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
int end(){
printf (" allways free memory (deallocate ) to avoid memory leaks \n"); // https://en.wikipedia.org/wiki/C_dynamic_memory_allocation
free (data);
free(edge);
PrintInfoAboutProgam();
return 0;
}
// ********************************************************************************************************************
/* ----------------------------------------- main -------------------------------------------------------------*/
// ********************************************************************************************************************
int main () {
setup ();
//DrawImagerOfBD(data);
//SaveArray2PGMFile (data, iWidth+0.2, "BD/J");
//ComputeBoundaries(data, edge);
//SaveArray2PGMFile (edge, iWidth+0.3, "boundaries of BD/J");
DrawImagerOfLSM(data);
SaveArray2PGMFile (data, iWidth+scale +0.6, "LSM/J");
ComputeBoundaries(data, edge);
SaveArray2PGMFile (edge, iWidth+scale+0.7, "boundaries of LSM/J");
CopyBoundaries(edge, data);
SaveArray2PGMFile (data, iWidth+scale+0.8, "LSM/J with boundaries");
end();
return 0;
}
- ↑ The Science of Fractal Images by Michael F. Barnsley, Robert L. Devaney, Benoit B. Mandelbrot, Heinz-Otto Peitgen, Dietmar Saupe, Richard F. Voss
Historia pliku
Kliknij na datę/czas, aby zobaczyć, jak plik wyglądał w tym czasie.
| Data i czas | Miniatura | Wymiary | Użytkownik | Opis | |
|---|---|---|---|---|---|
| aktualny | 19:24, 12 kwi 2020 | | 2000 × 2000 (646 KB) | wikimediacommons>Soul windsurfer | Uploaded own work with UploadWizard |
Lokalne wykorzystanie pliku
Poniższa strona korzysta z tego pliku:
