Here is my code. It shuld give similar Julia image in gray
For interior it gives similar result, but exterior not
How can I improve it?
/*
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 >i.txt
time ./a.out >e.txt
=======================
# gnuplot "i.plt"
set terminal svg enhanced background rgb 'white'
set xlabel "re(z)"
set ylabel "DLD"
set title "Relation between z and DLD in the interior of Julia set for c = -1"
set output "interior.svg"
plot "i.txt" with lines
----------------------
d0 - db = 5.0000000000000000 - 4.5389870050569598 = 0.4610129949430402
allways free memory (deallocate ) to avoid memory leaks
Numerical approximation of Julia set for fc(z)= z^2 + c
parameter c = ( -1.0000000000000000 ; 0.0000000000000000 )
Image Width = 4.000000 in world coordinate
PixelWidth = 0.004004
Maximal number of iterations = iterMax = 1000
ratio of image = 1.000000 ; it should be 1.000 ...
gcc version: 7.5.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
// https://sourceforge.net/p/predef/wiki/Standards/
#if defined(__STDC__)
# define PREDEF_STANDARD_C_1989
# if defined(__STDC_VERSION__)
# if (__STDC_VERSION__ >= 199409L)
# define PREDEF_STANDARD_C_1994
# endif
# if (__STDC_VERSION__ >= 199901L)
# define PREDEF_STANDARD_C_1999
# endif
# endif
#endif
/* --------------------------------- 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 = 5000; //
// 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.8; //-0.05;
static const double ZxMax = 1.8; //0.75;
static const double ZyMin = -1.8; //-0.1;
static const double ZyMax = 1.8; //0.7;
static double PixelWidth; // =(ZxMax-ZxMin)/ixMax;
static double PixelHeight; // =(ZyMax-ZyMin)/iyMax;
static double ratio;
// complex numbers of parametr plane
//https://fractalforums.org/code-snippets-fragments/74/lagrangian-descriptors-fragment-code/3612/msg22426#msg22426
double complex c = -1.05204872; // parameter of function fc(z)=z^2 + c
// attracting period 2 cycle
double complex z21 = 0.049822582293598;
double complex z22 = -1.049822582293598;
double complex z1a = -0.641073494565534; // alfa
complex double z1b = 1.641073494565534; // beta
double complex zb = 0.203208556149733 +0.377005347593583*I; // point on the boubndary of main componnet
int ARe = 60;
int ERe = 16;
double ER ;//= 1e60;
double AR ;//= 1e-16; // bigger values do not works
// DLD
const int N = 20; // fixed number : maximal number of iterations
double p = 0.01444322; //
// DLD colors
double me = 1.0;
double mi = 10.0;
double d21; // = lagrangian(z21, c, N, p);
double d22; // = lagrangian(z22, c, N, p);
double db;// = lagrangian(z1a, c, N, p);
double dd;// = d1a-d21;
/* colors = shades of gray from 0 to 255 */
unsigned char iColorOfExterior = 250;
unsigned char iColorOfInterior = 200;
unsigned char iColorOfBoundary = 0;
/* ------------------------------------------ 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) ****************************
/* ----------- 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;
}
// ***************************************************************************************************************************
// ************************** DLD/J*****************************************
// ****************************************************************************************************************************
/* partial pnorm
input: z , zn = f(z), p
output ppn
*/
double ppnorm( complex double z, complex double zn, double p){
double s[2][3]; // array for 2 points on the Riemann sphere
int j;
double d; // denominator
double x;
double y;
double ds;
double ppn = 0.0;
// map from complex plane to riemann sphere
// z
x = creal(z);
y = cimag(z);
d = x*x + y*y + 1.0;
s[0][0] = (2.0*x)/d;
s[0][1] = (2.0*y)/d;
s[0][2] = (d-2.0)/d; // (x^2 + y^2 - 1)/d
// zn
x = creal(zn);
y = cimag(zn);
d = x*x + y*y + 1.0;
s[1][0] = (2.0*x)/d;
s[1][1] = (2.0*y)/d;
s[1][2] = (d-2.0)/d; // (x^2 + y^2 - 1)/d
// sum
for (j=0; j <3; ++j){
ds = fabs(s[1][j] - s[0][j]);
// normal: neither zero, subnormal, infinite, nor NaN
//if (fpclassify (ds) !=FP_INFINITE)
//if (isnormal(ds))
// it is solved by if (cabs(z) > 1e60 ) break; procedure in parent function
ppn += pow(ds,p); // |ds|^p
// else {ppn = 10000.0; printf("ds = infty\t");} //
}
return ppn;
}
// DLD = Discret Lagrangian Descriptior
double lagrangian( complex double z0, complex double c, int iMax, double p ){
int i; // number of iteration
double d = 0.0; // DLD = sum
double ppn; // partial pnorm
complex double z = z0;
complex double zn; // next z
for (i=0; i<iMax; ++i){
zn = z*z +c; // complex iteration
ppn = ppnorm(z, zn, p);
d += ppn; // sum
//
z = zn;
//if (! isnormal(d)) { return 0.0; } // not works
if (cabs(z) > ER ) break; // exterior : big values produces NAN error in ppnorm computing
if (cabs(z - z21) < AR || cabs(z -z22)< AR)
{ // interior
//return FP_ZERO;
d = -d;
break;
}
}
if (d<0.0) {// interior
// d(z1a) - d(z21) = -0.0804163521959989
d = - d;
d = (db - d) /dd ; // normalize, see test_interior
//if (d>1.0) {printf("d int > 1.0\n");d = d - (int)d;} // fractional part}
}
else {
d = d/((double)i); // averaging not summation
d = d*me;} // exterior
return d;
}
unsigned char ComputeColorOfDLD(complex double z){
//double cabsz;
int iColor;
double d;
d = lagrangian(z,c, N,p);
//if ( d==FP_ZERO)
// {iColor = 255;}
// else
iColor = (int)(d*255) % 255; // nMax or lower walues in denominator
return (unsigned char) iColor;
}
// plots raster point (ix,iy)
int DrawPointOfDLD (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 = ComputeColorOfDLD(z);
A[i] = iColor ; // interior
return 0;
}
// fill array
// uses global var : ...
// scanning complex plane
int DrawImagerOfDLD (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)
DrawPointOfDLD(A, ix, iy); //
}
return 0;
}
// test how values changes to tune color
// from z21 to zb
int test_interior(){
complex double z = z21;
complex double dz = cabs(z21 - zb)/15.0;
double r = cabs(zb);
printf("# re(z) \t d\n"); // gnuplot
while (creal(z) < r){ // from z21 to z1a
double d = lagrangian(z, c, N, p);
//int iColor = ComputeColorOfDLD(z);
// printf(" z = %.16f d = %.16f color = %d \n",creal(z), d, iColor);
printf(" %.16f %.16f \n",creal(z), d); // gnuplot
z += dz;
}
//
d21 = lagrangian(z21, c, N, p);
db = lagrangian(zb, c, N, p);
d22 = lagrangian(z22, c, N, p);
dd = db - d21;
printf("d(z21) = %.16f \n",d21);
printf("d(z22) = %.16f \n",d22);
printf("d(zb) = %.16f \n",db);
printf("d(zb) - d(z21) = %.16f\n",dd);;
return 0;
}
// test how values changes to tune color
int test_exterior(){
complex double z;
complex double z0 = 1.62;
complex double z1 = 3.0;
complex double dz = 0.001;
z = z0;
printf("# z d\n"); // gnuplot
while (creal(z) < creal(z1)){
double d = lagrangian(z, c, N, p);
//int iColor = ComputeColorOfDLD(z);
// printf(" z = %.16f d = %.16f color = %d \n",creal(z), d, iColor);
printf(" %.16f %.16f \n",creal(z), d); // gnuplot
z += dz;
}
//
double d0 = lagrangian(z0, c, N, p);
double d1 = lagrangian(z1, c, N, p);
double dd = d0 - d1;
printf("d0 - d1 = %.16f - %.16f = %.16f\n",d0, d1, dd);
return 0;
}
// *******************************************************************************************
// ********************************** save A array to pgm file ****************************
// *********************************************************************************************
int SaveArray2PGMFile( unsigned char A[], int a, int b, int c, int d,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, "%d_%d_%d_%d", a,b,c,d); /* */
char *filename =strcat(name,".pgm");
char long_comment[200];
sprintf (long_comment, "%s\tER = %e\tAR =%e", comment, ER, AR);
// save image array 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", long_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 (long_comment == NULL || strlen(long_comment) ==0)
printf("\n");
else printf (". Comment = %s \n", long_comment);
return 0;
}
int PrintCInfo(){
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 displayed in the console : export OMP_DISPLAY_ENV="TRUE"
printf("__STDC__ = %d\n", __STDC__ );
printf("__STDC_VERSION__ = %ld\n", __STDC_VERSION__ );
printf("c dialect = ");
switch(__STDC_VERSION__) { // the format YYYYMM
case 199409L : printf("C94\n"); break;
case 199901L : printf("C99\n"); break;
case 201112L : printf("C11\n"); break;
case 201710L : printf("C18\n"); break;
//default : /* Optional */
}
return 0;
}
int PrintProgramInfo()
{
// display info messages
printf ("Numerical approximation of Julia set for fc(z)= z^2 + c \n");
//printf ("iPeriodParent = %d \n", iPeriodParent);
//printf ("iPeriodOfChild = %d \n", iPeriodChild);
printf ("parameter c = ( %.16f ; %.16f ) \n", creal(c), cimag(c));
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", N);
printf ("ratio of image = %f ; it should be 1.000 ...\n", ratio);
//
return 0;
}
// *****************************************************************************
//;;;;;;;;;;;;;;;;;;;;;; setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
// **************************************************************************************
int setup ()
{
printf ("setup start\n");
/* 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 ...
ER = pow(10,ARe);
AR = pow(10,-ERe);
/* create dynamic 1D arrays for colors ( shades of gray ) */
data = malloc (iSize * sizeof (unsigned char));
if (data == NULL ){
fprintf (stderr, " Could not allocate memory");
return 1;
}
test_interior();
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);
PrintProgramInfo();
PrintCInfo();
return 0;
}
// ********************************************************************************************************************
/* ----------------------------------------- main -------------------------------------------------------------*/
// ********************************************************************************************************************
int main () {
setup ();
DrawImagerOfDLD(data);
SaveArray2PGMFile (data, iWidth,N,ERe, ARe,"DLD/J , name = iWidth_N_ER_AR");
//test_exterior();
end();
return 0;
}