/*

  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 =-1.0; //   {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;
}