#include <Magick++.h>
#include <math.h>
#include <stdlib.h>
#include <fstream.h>
#include <complex.h>
#include <sys/time.h>
#include <unistd.h>


// always remeber , the no. of panels always remain the same in destination image ,
// only the panel size can vary


/* _____________________-
things to store in file :
1. image size - x , y ( int )
2. panelsize 

3. the Fractal codes for each destination panel
*/

using namespace std;
using namespace Magick;

const long MAXPANELS=200000;
const int nsforms=8;
const double _gamma=0.75;

struct panel
{
  int x,y;
  int panelsize;
  double mean;
};

struct fractalcode
{
  unsigned char  x,y;
  char s;
  unsigned char  gamma;
  short int  beta;
};

struct sform
{
  int a,b,c,d;
  double con;
  int xpixels[1025];
  int ypixels[1025];
 
};

char infile[80],outfile[80];
double image[1024][1024];
double rimage[1024][1024];
int rows,cols;
int panelsize=4; // default
int dxpanels,dypanels;
panel destpanels[MAXPANELS];

int rxpanels,rypanels;
panel refpanels[MAXPANELS];
int rrows,rcols;

fractalcode fractalcodes[MAXPANELS];
sform sformlist[nsforms]={{1,0,0,1,0.5},
			  {-1,0,0,1,0.5},
			  {1,0,0,-1,0.5},
			  {-1,0,0,-1,0.5},
			  {0,1,1,0,0.5},
			  {0,-1,1,0,0.5},
			  {0,1,-1,0,0.5},
			  {0,-1,-1,0,0.5}};




//___fn. to get pixel for a panel
inline double getpixel(panel p,int x,int y)
{
  return image[p.y+y][p.x+x];
}

//__fn to set pixel for a panel value
inline void setpixel(panel p,int x,int y,double value)
{
  image[p.y+y][p.x+x]=value;
}

void prepareall(int panelsize)
{

  int i,j,k;
  
  int x,y;
  for(k=0;k<8;k++)
    {
      
      for(x=0;x<panelsize;x++)
	for(y=0;y<panelsize;y++)
	  {
	    int index=x+y*panelsize;
	       
	    if(sformlist[k].a==1)
	      sformlist[k].xpixels[index]=x;
	    else if(sformlist[k].a==-1)
	      sformlist[k].xpixels[index]=panelsize-1-x;
	    
	    if(sformlist[k].b==1)
	      sformlist[k].xpixels[index] =y;
	    else if(sformlist[k].b==-1)
	      sformlist[k].xpixels[index]=panelsize-1-y;
	    
	    if(sformlist[k].c==1)
	      sformlist[k].ypixels[index]=x;
	    else if(sformlist[k].c==-1)
	      sformlist[k].ypixels[index] =panelsize-1-x;
	    
	    if(sformlist[k].d==1)
	      sformlist[k].ypixels[index]=y;
	    else if(sformlist[k].d==-1)
	      sformlist[k].ypixels[index]=panelsize-1-y;
	    
	  }
    }
  

  cout << "\nprepareall finished..." << flush;
  
}



//___fn to get pixels in transformed ref panels
// get transformed pixel
inline double getrpixel(panel p,int x,int y,int s)
{
  /*  int xc,yc;
  sform sf=sformlist[s];
  
  if(sf.a==1)
    xc=x;
  else if(sf.a==-1)
    xc=panelsize-1-x;
  
  if(sf.b==1)
    xc=y;
  else if(sf.b==-1)
    xc=panelsize-1-y;
   
  if(sf.c==1)
    yc=x;
  else if(sf.c==-1)
    yc=panelsize-1-x;
    
  if(sf.d==1)
    yc=y;
  else if(sf.d==-1)
    yc=panelsize-1-y;
  */

  int index=x+y*panelsize;
  int x1=sformlist[s].ypixels[index];
  int x2=sformlist[s].xpixels[index];
  return rimage[p.y+x1][p.x+x2];
}


double getmean(panel p)
{
  int i,j,k;
  double sum=0;
  for(i=p.y;i<p.y+panelsize;i++)
    for(j=p.x;j<p.x+panelsize;j++)
      {
	sum+=image[i][j];
      }

  return sum/(panelsize*panelsize);
  
}

double getrmean(panel p)
{
  int i,j,k;
  double sum=0;
  for(i=p.y;i<p.y+panelsize;i++)
    for(j=p.x;j<p.x+panelsize;j++)
      {
	sum+=rimage[i][j];
      }

  return sum/(panelsize*panelsize);
  
}



// prepare the destination panels of panelsize - original image
void preparedestpanels()
{

  cout << "\npreparing destination panels..." << flush;
  
  int i,j,k;
  int count=0;
  dxpanels=cols/panelsize;
  dypanels=rows/panelsize;
  int ndestpanels=dxpanels*dypanels;
  
  for(i=0;i<rows;i+=panelsize)
    for(j=0;j<cols;j+=panelsize)
      {
	panel r;
	r.x=j; r.y=i;
	r.panelsize=panelsize;
	r.mean=getmean(r);
	
	destpanels[count++]=r;	
      }

  cout << "\nprepatedestpanels : generated " << count << " panels..." << flush;
  
}


// prepare the reference image which is half the size  of original image
// also prepare the reference panels
void preparerefimage()
{

  cout << "\npreparing reference image ... " << flush;
  
  int i,j,k;
  rrows=rows/2;
  rcols=cols/2;

  for(i=0;i<rcols;i++)
    for(j=0;j<rrows;j++)
      {
	int xx=2*i;
	int yy=2*j;

	double val=(image[yy][xx]+image[yy+1][xx]+image[yy][xx+1]+image[yy+1][xx+1]);
	rimage[j][i]=(val)/4;
      }

  cout << "\nprepared the ref image by averaging the destination image..." << flush;
  

  rxpanels=rcols-(panelsize-1);
  rypanels=rrows-(panelsize-1);

  
  int nrefpanels=rxpanels*rypanels;

  cout << "\npreparerefimage : no of panels - " << nrefpanels << flush;
  
  // prepare the reference panels
  int count=0;
  for(i=0;i<rypanels;i++)
    for(j=0;j<rxpanels;j++)
      {
	
	panel r;
	r.x=j;r.y=i;
	r.panelsize=panelsize;
	r.mean=getrmean(r); // get mean
	refpanels[count++]=r;	
      }  

  cout << "\npreparerefimage : prepared ref panels. count : " << count << flush;
}


void readimage(char* filename)
{


  cout << "\nreadimage : " << filename << flush;
  Image source;
  
  source = Image(filename);
  cols = source.baseColumns();  
  rows = source.baseRows();
  source.classType(DirectClass);
  source.quantizeColorSpace(GRAYColorspace);

  Pixels viewsrc(source);
  PixelPacket *pixelsrc = viewsrc.get(0,0,cols,rows);
  
  for(int x=0;x<rows;x++)
    for(int y=0;y<cols;y++)
      {
	image[x][y] = ((ColorGray)(*pixelsrc)).shade();	
	pixelsrc++;
      }
  pixelsrc-= cols*rows;
  
}

void writeimage(char* filename,int ask=0)
{

  cout << "\nwriteimage : " << filename << flush;
  
  char buffer[80];
  if(ask)
    {

      char choice;
      cout << "\nWant to save the image(y/n)  : ";
      cin >> choice;
      if(choice!='y' && choice!='Y')
	return;

      cout << "Enter filename : " ;
      cin >> buffer;
    }
  
  Image fspectrum(Geometry(cols,rows),"white");
  fspectrum.classType(DirectClass);
  fspectrum.quantizeColorSpace(GRAYColorspace);
  Pixels fview(fspectrum);
  PixelPacket *fpixels = fview.get(0,0,cols,rows);

  for(int x=0;x<rows;x++)
    for(int y=0;y<cols;y++)
      {
	(*fpixels) = ColorGray(image[x][y]);
	fpixels++;
      }
    
  fview.sync();
  fspectrum.write((ask==1) ? buffer : filename);

}

//___parameters are panelsize only now
void compress(char* filename)
{

  cout << "\nCompressing image with panelsize : " << panelsize << flush;
  // panelsize=4;

  prepareall(panelsize); // prepare the s form transformation pixel array
  
  preparedestpanels();
  preparerefimage();
  
  int ndestpanels=dxpanels*dypanels;
  int nrefpanels=rxpanels*rypanels;
  
  double error,beta,gamma,diff;
  double besterror,bestbeta,bestgamma,bestform;
  int bestsform;
  panel bestrpanel;
  
  int i,j,k,x,y;
  cout << "\ncompress:  no of destination panels :: " << ndestpanels << flush;

  double dr,d,r,d2,t;
  for(i=0;i<ndestpanels;i++) // loop over all destpanels
    {
      besterror=9999999;

      r=0;

      // compute statistics
      
      for(x=0;x<panelsize;x++)
	for(y=0;y<panelsize;y++)
	  r+=image[destpanels[i].y+y][destpanels[i].x+x];

      
      for(j=0;j<nrefpanels;j++) // search for the best refpanel
	{
	  // calculate difference in mean;
	  //beta=destpanels[i].mean-refpanels[j].mean;


	  // compute alpha and beta

	  // try all orientations of ref image to suit best image for this destination panel
	  int s;
	  
	   for(s=0;s<nsforms;s++)
	    {
	      error=0;
	      

	      dr=0;
	      d=0;
	      d2=0;
	      // calculate metrics
	      for(x=0;x<panelsize;x++)
		for(y=0;y<panelsize;y++)
		  {
		    t=getrpixel(refpanels[j],x,y,s);
		    d+=t;
		    d2+=t*t;
		    dr+=image[destpanels[i].y+y][destpanels[i].x+x]*getrpixel(refpanels[j],x,y,s);

		  }

	      //r*=255;
	      //d*=255;
	      //d2*=(255*255);
	      //dr*=(255*255);
	      
	      gamma=(panelsize*panelsize*dr-d*r)/(panelsize*panelsize*d2-d*d);
	      if(gamma>1.0) gamma=1.0;
	      if(gamma<-1.0) gamma=-1.0;
	      
	      //beta=(d*d-r*r)/(panelsize*panelsize*d2-d*d);

	      //cout << "\nGamma for refpanel : " << j << " = " << (panelsize*panelsize*dr-d*r) << " / " << (panelsize*panelsize*d2-d*d) << "=" << gamma;

	      // calculate beta

	      beta=0;
	      for(x=0;x<panelsize;x++)
		for(y=0;y<panelsize;y++)
		  beta+=image[destpanels[i].y+y][destpanels[i].x+x]-gamma*getrpixel(refpanels[j],x,y,s);
	      
	      
	      beta/=(panelsize*panelsize);
								     
		 
	      for(x=0;x<panelsize;x++)
		for(y=0;y<panelsize;y++)
		  {
		    
		    // get diff between corresponding pixels
		    diff=image[destpanels[i].y+y][destpanels[i].x+x]-(gamma*getrpixel(refpanels[j],x,y,s)+beta);
		    // calc variance 
		    error+=(diff)*(diff)*255*255;
		  }
	  
		  
	      if(error<besterror)
		{
		  besterror=error;
		  bestsform=s;
		  bestbeta=beta;
		  bestgamma=gamma;
		  bestrpanel=refpanels[j];
		}
	    } // for s
	  

	   if(besterror==0) break;
	} // for refpanels

      // add the best region panel for dest panel here
      
      fractalcodes[i].x=(unsigned char)bestrpanel.x;
      fractalcodes[i].y=(unsigned char)bestrpanel.y;
      fractalcodes[i].s=bestsform;
      fractalcodes[i].beta=(bestbeta*255);
      fractalcodes[i].gamma=(unsigned char)(bestgamma*127+128.0);
      cout << "\ndone for destination panel : " << i << " out of "  << ndestpanels  << flush;
      //cout << "\nDestpanels coords : " << destpanels[i].x << "," << destpanels[i].y << flush;
      //cout << "\nbestrefpanel coords : " << bestrpanel.x << "," << bestrpanel.y << " Best gamma : " << bestgamma<<"-"<<int(fractalcodes[i].gamma-128)/127.0 << " best beta : " << bestbeta <<"-"<< int(fractalcodes[i].beta)/255.0 <<  flush;
      //cout << "\n  Error : " << error << flush;
      
    }
  

  // compression done
  // now fractal codes have the fractal information
  // save to file
  
  cout << "\ncompress : writing the fractal codes to filename : " << filename;

  FILE* fp=fopen(filename,"wb");
  fwrite(&cols,sizeof(cols),1,fp);
  fwrite(&rows,sizeof(rows),1,fp);

  //  fwrite(&dxpanels,sizeof(dxpanels),1,fp);
  fwrite(&panelsize,sizeof(panelsize),1,fp);

  fwrite(fractalcodes,sizeof(fractalcode),ndestpanels,fp);

  fclose(fp);

  cout  << "\nCompress : All the records written to file..." << flush;
    
}



// __xsize , ysize , xpanels, ypanels --- from original image parameters
void decompress(char* filename,int xsize,int ysize,int psize)
{ 
  cout << "\ndecompress called..." << flush;
  
  //___routine for reading parameters from file____

  int xpanels,ypanels;
    
  xpanels=xsize/psize;
  ypanels=ysize/psize; // always same


  int ncodes=xpanels*ypanels;

  //  cout << "\ndecompress : no of panels in destination image : " << ncodes << flush;
  
  //  cout  << "\ndecompress:  read the parameters from file ... now decompressing the image ..." << flush;
  
  double scalefactor=cols/xsize;

  int ndestpanels=xpanels*ypanels; // always the same
  
  
  panelsize=cols/xpanels;
  
  prepareall(panelsize);
  
  preparedestpanels();
  preparerefimage();
  
 
  int i,j,k;

  cout << "\nScale factor : " << scalefactor << "\npanelsize : "<<panelsize << "\n" << flush;
  
  for(i=0;i<ndestpanels;i++)
    {

      //cout << "\ndecompress : trying to decompress panel no. : " << i << flush;
      
      int x=(int)(scalefactor*int(fractalcodes[i].x)+0.49);
      int y=(int)(scalefactor*int(fractalcodes[i].y)+0.49);
      if(x>=rxpanels) {x=rxpanels-1;cout << "x>=rxpanels..." << flush;}
      if(y>=rypanels) {y=rypanels-1;cout << "y>=rypanels..." << flush;}

      //cout << "\nRefpanel no  : " << x+y*rxpanels << " sform : " << int(fractalcodes[i].s) << "\ndestpanel no : " << i << flush;
      panel refpanel=refpanels[x+y*rxpanels];
      panel destpanel=destpanels[i];

      //cout << "\nStarting the double loop for transforming the doman to range block " << flush;
      for(j=0;j<panelsize;j++)
	for(k=0;k<panelsize;k++)
	  {
	    //	    cout << "\ncalling getrpixel... value of s : " << fractalcodes[i].s << flush;
	    double val=getrpixel(refpanel,j,k,int(fractalcodes[i].s))*double(int(fractalcodes[i].gamma)-128)/127.0+double(int(fractalcodes[i].beta))/255.0;
	    //double val=rimage[refpanel.y+k][refpanel.x+j]+double(fractalcodes[i].beta)/255.0;
	    
	    if(val<=0)
	      {val=0;
	      //cout << "\nval<0..." << flush;
	      }
	    else if
	      (val>=1) {val=1;
	      //cout << "\nval>1" << flush;
	    }

	    // cout << "\ncalling setpixel at ..." << j << "," <<  k  << " Destpanel : " << destpanel.x << "," << destpanel.y << flush;
	    setpixel(destpanel,j,k,val);    
	  }
    }
  
  // one iteration is finished here
  // repeat this many times to get the final image
  }

void main()
{
  while(1)
    {
      

      cout << endl << "    MAIN MENU ";
      cout << endl << "---------------";
      cout << endl << "1. Compress image";
      cout << endl << "2. Decompress image";

      cout << endl << "3. Exit";
  
      cout << "\n\n\n";

      cout << "Enter choice :";

      char buffer[80];
      fgets(buffer,79,stdin);
      int choice=atoi(buffer);

      if(choice==3)
	break;
      if(choice<1 || choice>>3)
	continue;

 
      cout  << "\n\nEnter infile : ";
      cin >> infile;
      cout << "\nEnter outfile :";
      cin >> outfile;

      if(choice==1) // compress file
	{
	  cout << "\nPlease enter the panel size : ";
	  cin >> panelsize;
	  
	  readimage(infile);
	  if(cols%panelsize!=0 || rows%panelsize!=0)
	    {
	      cout << "\nImage resolution is not a multiple of panelsize. Re-enter." << flush;
	      continue;
	    }
	  
	  
	  compress(outfile);
	  
	  cout << "\nMain : succesfully returned from function compress " << flush;
	}
      
      if(choice==2) // decompress image
	{
	  int i,j;
	  
	  // first read the square image
	  //int x,y;
	  //cout << "\n Enter xres , yres :" << flush;
	  //cin >> x,y;

	  int xsize,ysize,psize;
	  FILE* fp=fopen(infile,"rb");
	  fread(&xsize,sizeof(xsize),1,fp);
	  fread(&ysize,sizeof(ysize),1,fp);
	  fread(&psize,sizeof(psize),1,fp);

	  int xpanels,ypanels;
	  
	  xpanels=xsize/psize;
	  ypanels=ysize/psize; // always same
	  int ncodes=xpanels*ypanels;

	  
	  fread(fractalcodes,sizeof(fractalcode),ncodes,fp);
	  
	  fclose(fp);

	  cout << "\nCompression Panel size is " << psize << endl;
	  cout << "\nEnter output panelsize : "  << flush;
	  cin >> panelsize;
	  

	  int x,y;
	  x=(xsize/psize)*panelsize;
	  y=(ysize/psize)*panelsize;

	  
	  cout << "\nInput resolution : " << xsize << "," << ysize ;
	  cout << "\nOutput resolution : " << x << "," << y;
	  
	  /*
	  char buffer[80];
	  sprintf(buffer,"convert -scale %dx%d square.gif s.gif",x,y);
	  system(buffer);
	  */

	  
	  cols=x;rows=y;
	  for(i=0;i<rows;i++)
	    for(j=0;j<cols;j++)
	      if(j>cols/2) image[i][j]=1; else image[i][j]=0;

	  writeimage("s.gif");

	  
	  cout << "\nReading s.gif ( the square image )..." << flush;
	  readimage("s.gif");
	  // this has gone to image array
	  
	  
	  // 8 iterations of decompress
	  for(i=0;i<10;i++)
	    {
	      cout << "\nCalling iteration no : "<< i << flush;
	      decompress(infile,xsize,ysize,psize);
	      sprintf(buffer,"file%d.gif",i);
	      writeimage(buffer);
	      
	    }
	  
	  // now image is in image
	  
	  cout << "\nWriting the image to : " << outfile << flush;
	  
	  writeimage(outfile);
	  cout << "\nWritten the image to  : " << outfile << flush;
	  
	  cout << "\nCheck it out..." << flush;
	  
	  
	}
      
    }
  
}