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


using namespace Magick;

class CImage
{
public:
  int maxX;
  int maxY;
  int cRes;
  int **intensity;

  int loadImage(char*);
  int saveImage(char*);
  void saveHistogram(char*);
  
  void negative();
  void contrastStretch(double r1, double r2, double s1, double s2);
  void averaging(int mask);
  void median(int mask);
  void localHistEq(int mask);
  void histSpecification(char*);
  void highPassFilter();
  void highBoost(int mask,float A, bool average);
  void globalHistEq();
  void print();
  void showHistogram();
  void showImage();
  void MEqualize(char* file);
  void MMedianFilter(char* file, int mask);
  
  CImage()
  {
    maxX = maxY = 0;
    intensity = NULL;
    cRes = 255;
  }

  CImage(int maxx, int maxy)
  {
    maxX = maxx;
    maxY = maxy;
    intensity = new int*[maxY];
    for(int i=0;i<maxY;i++)
      intensity[i] = new int[maxX];
    cRes = 255;
  }
  ~CImage()
  {
    int i;
    for(i=0;i<maxY;i++)
      delete intensity[i];
    delete [] intensity;
  }
};

void CImage::MEqualize(char* file)
{
  struct timeval tv1,tv2;
  struct timezone tz1,tz2;
  long timegap;
  Image image(file);
  gettimeofday(&tv1,&tz1);
  image.equalize();
  gettimeofday(&tv2,&tz2);
  timegap = (tv2.tv_sec - tv1.tv_sec) * 1000000 + tv2.tv_usec - tv1.tv_usec;
  cout<<"Time taken : "<<timegap<<" microseconds.\n";
  image.display();
  image.write("temp1.jpg");
}

void CImage::MMedianFilter(char* file, int mask)
{
  struct timeval tv1,tv2;
  struct timezone tz1,tz2;
  long timegap;
  Image image(file);
  gettimeofday(&tv1,&tz1);
  image.medianFilter(mask/2);
  gettimeofday(&tv2,&tz2);
  timegap = (tv2.tv_sec - tv1.tv_sec) * 1000000 + tv2.tv_usec - tv1.tv_usec;
  cout<<"Time taken : "<<timegap<<" microseconds.\n";
  image.display();
  image.write("temp2.jpg");
}

void CImage::negative()
{
  int i,j;
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      intensity[i][j] = cRes - intensity[i][j];
}

int CImage::loadImage(char* file)
{
  Image image(file);
  maxX=image.columns();
  maxY=image.rows();
  image.classType(DirectClass);
  Pixels view(image);
  PixelPacket *pixels=view.get(0,0,maxX,maxY);
  int i,j;
  intensity = new int*[maxY];
  for(int i=0;i<maxY;i++)
    intensity[i] = new int[maxX];
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	cout.flush();
        ColorGray gray=Color(*(pixels+i*maxX+j));
        intensity[i][j]=(int)(gray.shade()*cRes);
      }
}

int CImage::saveImage(char* file)
{
  Image image(Geometry(maxX,maxY), "white");
  image.classType(DirectClass);
  Pixels view(image);
  PixelPacket *pixels=view.get(0,0,maxX,maxY);
  int i,j;
  for(int i=0;i<maxY;i++)
    for(int j=0;j<maxX;j++)
      (*(pixels+i*maxX+j))=ColorGray(float(intensity[i][j])/(float)cRes);
  image.write(file);
}

void CImage::print()
{
  for(int i=0;i<maxY;i++)
    {
      for(int j=0;j<maxX;j++)
	cout<<" "<<intensity[i][j];
      cout<<"\n";
    }
}

void CImage::contrastStretch(double r1, double r2, double s1, double s2)
{
  for(int i=0;i<maxY;i++)
    for(int j=0;j<maxX;j++)
      {
	if(intensity[i][j]<r1*cRes)
	  intensity[i][j] = (int)(((float)intensity[i][j])*s1/r1);
	else if(intensity[i][j]<r2*cRes)
	  intensity[i][j] = (int)(s1*cRes + (s2-s1)/(r2-r1)*float(intensity[i][j]-r1*cRes));
	else
	  intensity[i][j] = (int)(s2*cRes + (float)(1.0-s2)/float(1.0-r2)*float(intensity[i][j]-r2*cRes));
      }
}

void CImage::highPassFilter()
{
  int mask = 3;
  int** newIntensity = new int*[maxY];
  for(int i=0;i<maxY;i++)
    newIntensity[i] = new int[maxX];
  int i,j;
  int back = -(mask/2);
  int forward = mask/2;
  int i1,j1;
  int totalVal,totalInt;
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	//cout<<i<<","<<j<<" ";cout.flush();
	totalVal = 0;
	totalInt = 0;
	for (i1 = i+back; i1 <= i+forward ; i1++)
          for (j1 = j+back; j1 <= j+forward ; j1++)
	    {
	      if((i1>=0)&&(j1>=0)&&(i1<maxY)&&(j1<maxX))
		{
		  totalInt -= intensity[i1][j1];
		  totalVal++;
		}
	    }
	totalInt += 9*intensity[i][j];
	//cout<<"a"<<intensity[i][j];cout.flush();	
	newIntensity[i][j] = totalInt/totalVal;
	//cout<<"b";cout.flush();
      }
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      intensity[i][j] = newIntensity[i][j];
  for(i=0;i<maxY;i++)
    delete newIntensity[i];
  delete [] newIntensity;
}

void CImage::highBoost(int mask,float A,bool average=true)
{
  int i,j;
  int** newIntensity = new int*[maxY];
  for(i=0;i<maxY;i++)
    newIntensity[i] = new int[maxX];
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      newIntensity[i][j] = intensity[i][j];
  if(average)
    averaging(mask);
  else
    median(mask);
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      intensity[i][j] = int(A * (float)newIntensity[i][j]) - intensity[i][j];  
  for(i=0;i<maxY;i++)
    delete newIntensity[i];
  delete [] newIntensity;
}

void CImage::averaging(int mask)
{
  int** newIntensity = new int*[maxY];
  for(int i=0;i<maxY;i++)
    newIntensity[i] = new int[maxX];
  int i,j;
  int back = -(mask/2);
  int forward = mask/2;
  int i1,j1;
  int totalVal,totalInt;
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	//cout<<i<<","<<j<<" ";cout.flush();
	totalVal = 0;
	totalInt = 0;
	for (i1 = i+back; i1 <= i+forward ; i1++)
          for (j1 = j+back; j1 <= j+forward ; j1++)
	    {
	      if((i1>=0)&&(j1>=0)&&(i1<maxY)&&(j1<maxX))
		{
		  totalInt += intensity[i1][j1];
		  totalVal++;
		}
	    }
	//cout<<"a"<<intensity[i][j];cout.flush();	
	newIntensity[i][j] = totalInt/totalVal;
	//cout<<"b";cout.flush();
      }
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      intensity[i][j] = newIntensity[i][j];
  for(i=0;i<maxY;i++)
    delete newIntensity[i];
  delete [] newIntensity;
}

void CImage::histSpecification(char* file)
{
  CImage c;
  c.loadImage(file);
  globalHistEq();
  int *val;
  val = new int[cRes+1];
  int total = 0;
  int i,j;
  
  for(i=0;i<=cRes;i++)
    val[i] = 0;
  for(i=0;i<c.maxY;i++)
    for(j=0;j<c.maxX;j++)
      {
	val[c.intensity[i][j]]++;
	total++;
      }

  int *newI;
  newI = new int[cRes+1];
  int *newJ;
  newJ = new int[cRes+1];
  int v = 0;
  for(i=0;i<=cRes;i++)
    {
      v += val[i];
      newI[i] = (int)((float)v/(float)total*cRes); 
    }

  i=0;
  int oldi= 0;
  int k = 0;
  while((i<cRes)&&(newI[i]==newI[i+1]))
    i++;
  while(i<=cRes)
    {
      for(j=k;j<=newI[i];j++)
	newJ[j] = oldi;
      k = newI[i]+1;
      i++;
      oldi = i;
      while((i<cRes)&&(newI[i]==newI[i+1]))
	i++;
  }
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      intensity[i][j] = newJ[intensity[i][j]];

  delete [] val; 

}


int compare(const void* a, const void* b)
{
  return(*((int*)a) > *((int*)b));
}

void CImage::median(int mask)
{
  int** newIntensity = new int*[maxY];
  for(int i=0;i<maxY;i++)
    newIntensity[i] = new int[maxX];
  int* val =  new int[mask*mask];
  int i,j;
  int back = -(mask/2);
  int forward = mask/2;
  int i1,j1;
  int totalVal;
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	//cout<<i<<","<<j<<" ";cout.flush();
	totalVal = 0;
	for (i1 = i+back; i1 <= i+forward ; i1++)
          for (j1 = j+back; j1 <= j+forward ; j1++)
	    {
	      if((i1>=0)&&(j1>=0)&&(i1<maxY)&&(j1<maxX))
		{
		  val[totalVal] = intensity[i1][j1];
		  totalVal++;
		}
	    }
	qsort(val,totalVal,sizeof(int),compare);
	
	//cout<<"a"<<intensity[i][j];cout.flush();	
	newIntensity[i][j] = val[totalVal/2];
	//cout<<"b";cout.flush();
      }
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      intensity[i][j] = newIntensity[i][j];
  for(i=0;i<maxY;i++)
    delete newIntensity[i];
  delete [] newIntensity;
}


void CImage::localHistEq(int mask)
{
  int** newIntensity = new int*[maxY];
  for(int i=0;i<maxY;i++)
    newIntensity[i] = new int[maxX];

  int *val;
  val = new int[cRes+1];
  int total = 0;
  int i,j;
  int back = -(mask/2);
  int forward = mask/2;
  int *newI;
  newI = new int[cRes+1];
  int v = 0;
  int i1,j1,i2,j2;
  for(i=0;i<maxY;i++)
    {
      total = 0;
      for(i2=0;i2<=cRes;i2++)
	val[i2] = 0;
      for (i1 = i+back; i1 <= i+forward ; i1++)
	{
	  for (j1 = back; j1 <= forward ; j1++)
	    {
	      if((i1>=0)&&(j1>=0)&&(i1<maxY)&&(j1<maxX))
		{
		  val[intensity[i1][j1]]++;
		  total++;
		}
	    }
	}
      v=0;
      for(i1=0;i1<=intensity[i][0];i1++)
	v += val[i1];
      newIntensity[i][0] = int(float(v)/float(total)*float(cRes));
      for(j=1;j<maxX;j++)
      {
	j2 = j + back - 1;
	for (i1 = i+back; i1 <= i+forward ; i1++)
	  {
	    if((i1>=0)&&(j2>=0)&&(i1<maxY)&&(j2<maxX))
	      {
		val[intensity[i1][j2]]--;
		total--;
	      }
	  }
	j2 = j+forward;
	for (i1 = i+back; i1 <= i+forward ; i1++)
	  {
	    if((i1>=0)&&(j2>=0)&&(i1<maxY)&&(j2<maxX))
	      {
		val[intensity[i1][j2]]++;
		total++;
	      }
	  }
	v=0;
	for(i1=0;i1<=intensity[i][j];i1++)
	  v += val[i1];
	newIntensity[i][j] = int(float(v)/float(total)*float(cRes));
	//cout<<"b";cout.flush();
      }
    }
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      intensity[i][j] = newIntensity[i][j];
  
  delete [] val;
}

void CImage::showImage()
{
  saveImage("temp.jpg");
  system("xv temp.jpg &");
}

void CImage::showHistogram()
{
  int *val;
  val = new int[cRes+1];
  int total = 0;
  int i,j;
  
  for(i=0;i<=cRes;i++)
    val[i] = 0;
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	val[intensity[i][j]]++;
	total++;
      }

  fstream fout;
  fout.open("data",ios::out);
  for(i=0;i<=cRes;i++)
    fout<<i<<" "<<val[i]<<"\n";
  fout.close();
  
  fout.open("plotfile",ios::out);
  fout<<"set size 0.5,0.5\n";
  fout<<"set xlabel \"Graylevel\"\nset ylabel \"Frequency\"\n";
  fout<<"plot[0:255] \"data\"  title \"Histogram\" with boxes\n";
  fout.close();

  system("gnuplot -persist plotfile &");

  delete [] val;
}

void CImage::saveHistogram(char* file)
{
  int *val;
  val = new int[cRes+1];
  int total = 0;
  int i,j;
  
  for(i=0;i<=cRes;i++)
    val[i] = 0;
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	val[intensity[i][j]]++;
	total++;
      }

  fstream fout;
  fout.open("data",ios::out);
  for(i=0;i<=cRes;i++)
    fout<<i<<" "<<val[i]<<"\n";
  fout.close();
  
  fout.open("plotfile",ios::out);
  fout<<"set size 0.5,0.5\n";
  fout<<"set terminal gif\n";
  fout<<"set output \""<<file<<"\"\n";
  fout<<"set xlabel \"Graylevel\"\nset ylabel \"Frequency\"\n";
  fout<<"plot[0:255] \"data\"  title \"Histogram\" with boxes\n";
  fout<<"plot \"data\" with boxes";
  fout.close();

  system("gnuplot plotfile &");

  delete [] val;
}

void CImage::globalHistEq()
{
  int *val;
  val = new int[cRes+1];
  int total = 0;
  int i,j;
  
  for(i=0;i<=cRes;i++)
    val[i] = 0;
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	val[intensity[i][j]]++;
	total++;
      }

  fstream fout;
  fout.open("beforeGHE",ios::out);
  for(i=0;i<=cRes;i++)
    fout<<i<<" "<<val[i]<<"\n";
  fout.close();

  int *newI;
  newI = new int[cRes+1];
  int v = 0;
  for(i=0;i<=cRes;i++)
    {
      v += val[i];
      newI[i] = (int)((float)v/(float)total*cRes); 
    }

  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      intensity[i][j] = newI[intensity[i][j]];

  for(i=0;i<=cRes;i++)
    val[i] = 0;
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      val[intensity[i][j]]++;

  fout.open("afterGHE",ios::out);
  for(i=0;i<=cRes;i++)
    fout<<i<<" "<<val[i]<<"\n";
  fout.close();
  delete [] val;
}

void main()
{
  CImage i;
  char fname[20];
  int option = 1;
  float r1,r2,s1,s2;
  int mask;
  float A;
  struct timeval tv1,tv2;
  struct timezone tz1,tz2;
  long timegap;

  while(option !=0)
    {
      cout<<"\n\n\n\n\tIMAGE TRANSFORMATION - SPACIAL\n\n";
      cout<<"1.  Load Image\n";
      cout<<"2.  View Image\n";
      cout<<"3.  Save Image\n";
      cout<<"4.  View Hisogram\n";
      cout<<"5.  Save Hisogram\n";
      cout<<"6.  Make Negative\n";
      cout<<"7.  Perform Contrast Streatching\n";
      cout<<"8.  Perform Global Histogram Equalization\n";
      cout<<"9.  Perform Local Histogram Equalization\n";
      cout<<"10. Perform Histogram Specification\n";
      cout<<"11. Perform Smoothing - Averaging Filter\n";
      cout<<"12. Perform Smoothing - Median Filter\n";
      cout<<"13. Perform Sharpening - High Boost Fiter\n";
      cout<<"14. Image Magick - Equalize\n";
      cout<<"15. Image Magick - MedianFilter\n";
      cout<<"0.  Exit\n\n";
      cout<<"Enter Option : ";
      cin>>option;
      
      switch(option)
	{
	case 1 :
	  cout<<"\n\nEnter name of image file (JPEG file) : ";
	  cin>>fname;
	  i.loadImage(fname);
	  cout<<"\n\nImage is loaded.\nSize is : "<<i.maxY<<"X"<<i.maxX<<"\n";
	  break;

	case 2 :
	  i.showImage();
	  break;

	case 3 :
	  cout<<"\n\nEnter name of image file (JPEG file) : ";
	  cin>>fname;
	  i.saveImage(fname);
	  cout<<"\n\nImage is saved.\n";
	  break;

	case 4 :
	  i.showHistogram();
	  break;

	case 5 :
	  cout<<"\n\nEnter name of image file (GIF file) : ";
	  cin>>fname;
	  i.saveHistogram(fname);
	  cout<<"\n\nImage is saved.\n";
	  break;

	case 6 :
	  i.negative();
	  cout<<"Negation is done";
	  break;
	  
	case 7 :
	  cout<<"Enter r1 : ";
	  cin>> r1;
	  cout<<"Enter r2 : ";
	  cin>> r2;
	  cout<<"Enter s1 : ";
	  cin>> s1;
	  cout<<"Enter s2 : ";
	  cin>> s2;
	  i.contrastStretch(r1,r2,s1,s2);
	  cout<<"Contrast Steaching is done";
	  break;

	case 8 :
	  gettimeofday(&tv1,&tz1);
	  i.globalHistEq();
	  gettimeofday(&tv2,&tz2);
          timegap = (tv2.tv_sec - tv1.tv_sec) * 1000000 + tv2.tv_usec - tv1.tv_usec;
          cout<<"Time taken : "<<timegap<<" microseconds.\n";
	  cout<<"Global Histogram Equalization is done";
	  break;

	case 9 :
	  cout<<"\n\n Enter mask : ";
	  cin>>mask;
          gettimeofday(&tv1,&tz1);
	  i.localHistEq(mask);
	  gettimeofday(&tv2,&tz2);
          timegap = (tv2.tv_sec - tv1.tv_sec) * 1000000 + tv2.tv_usec - tv1.tv_usec;
          cout<<"Time taken : "<<timegap<<" microseconds.\n";
	  cout<<"Local Histogram Equalization is done";
	  break;

	case 10 :
	  cout<<"\n\nEnter name of image file (JPEG file) : ";
	  cin>>fname;
	  i.histSpecification(fname);
	  break;
	  
	case 11 :
	  cout<<"\n\n Enter mask : ";
	  cin>>mask;
	  i.averaging(mask);
	  cout<<"Averaging is done";
	  break;
  
	case 12 :
	  cout<<"\n\n Enter mask : ";
	  cin>>mask;
	  gettimeofday(&tv1,&tz1);
	  i.median(mask);
	  gettimeofday(&tv2,&tz2);
          timegap = (tv2.tv_sec - tv1.tv_sec) * 1000000 + tv2.tv_usec - tv1.tv_usec;
          cout<<"Time taken : "<<timegap<<" microseconds.\n";
	  cout<<"Median Filtering is done";
	  break;

	case 13 :
	  cout<<"\n\n Enter mask : ";
	  cin>>mask;
	  cout<<"Enter A : ";
	  cin>>A;
	
	  i.highBoost(mask,A);
	  cout<<"Sharpening is done";
	  break;

	case 14 :
	  cout<<"\n\nEnter name of image file (JPEG file) : ";
	  cin>>fname;
	  i.MEqualize(fname);
	  cout<<"\n\nImage is Equalized\n";
	  break;

	case 15 :
	  cout<<"\n\nEnter name of image file (JPEG file) : ";
	  cin>>fname;
	  cout<<"\n\n Enter mask : ";
	  cin>>mask;
	  i.MMedianFilter(fname,mask);
	  cout<<"\n\nImage is Smoothened\n";
	  break;

	case 0 :
	  cout<<"\n\n\nProgram made by : \n Ashish Gupta (98130) & \n Ashish Gupta (98131) \n\n\n";
	  break;
	  
	default:
	  cout<<"\n\nWrong Option.";
	}
    }
}