#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 print();
  void showImage();
  void DetectB();
  void DetectBS();
  void LoG();
  void HoughTransformLine(int nt,int nc,int n1);
  void HoughTransformCircle(int nc,int n1);
  //  void HoughTransform(int nt,int nc);
  
  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::LoG()
{
  int maskL[3][3] = {0,1,0,1,-4,1,0,1,0};
  // int maskG[5][5] = {1,4,7,4,1,4,16,26,16,4,7,26,41,26,7,4,16,26,16,4,1,4,7,4,1};
   int maskG[5][5] = {2,4,5,4,2,4,9,12,9,4,5,12,15,12,5,4,9,12,9,4,2,4,5,4,2};
  int** newIntensity = new int*[maxY];
  int i,j,i1,j1;
  for( i=0;i<maxY;i++)
    newIntensity[i] = new int[maxX];
  int totalVal,totalInt;
  //Laplacian
  
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	totalVal = 0;
	totalInt = 0;
	for (i1 = i-1; i1 <= i+1 ; i1++)
          for (j1 = j-1; j1 <= j+1 ; j1++)
	    {
	      if((i1>=0)&&(j1>=0)&&(i1<maxY)&&(j1<maxX))
		{
		  totalInt += intensity[i1][j1]*maskL[i1-i+1][j1-j+1];
		  totalVal++;
		}
	    }
	newIntensity[i][j] = totalInt/totalVal;
      }

  // Gausian
  
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	totalVal = 0;
	totalInt = 0;
	for (i1 = i-2; i1 <= i+2 ; i1++)
          for (j1 = j-2; j1 <= j+2 ; j1++)
	    {
	      if((i1>=0)&&(j1>=0)&&(i1<maxY)&&(j1<maxX))
		{
		  totalInt += newIntensity[i1][j1]*maskG[i1-i+2][j1-j+2];
		  totalVal++;
		}
	    }
	intensity[i][j] = totalInt/totalVal;
      }
    DetectBS();
    for( i=0;i<maxY;i++)
      delete[] newIntensity[i];
    delete[] newIntensity;
}

void CImage::DetectBS()
{
  
  for(int i=0;i<maxY-1;i++)
    for(int j=0;j<maxX-1;j++)
      if(((intensity[i][j] >=0) &&(intensity[i+1][j] <0)) ||((intensity[i][j] >=0) &&(intensity[i][j+1] <0)))
	intensity[i][j] = 255;
      else
	intensity[i][j] = 0;
  for(int i=0;i<maxY;i++)
    {
      intensity[i][0] = intensity[i][1] = intensity[i][maxX-2] = intensity[i][maxX-1] = 0;
      intensity[i][2] = intensity[i][3] = intensity[i][maxX-3] = intensity[i][maxX-4] = 0;

    }
  for(int j=0;j<maxX;j++)
    {
      intensity[0][j] = intensity[1][j] = intensity[maxY-2][j] = intensity[maxY-1][j] = 0;
      intensity[2][j] = intensity[3][j] = intensity[maxY-3][j] = intensity[maxY-4][j] = 0;
    }
}
void CImage::DetectB()
{
  
  for(int i=0;i<maxY-1;i++)
    for(int j=0;j<maxX-1;j++)
      if( ((intensity[i][j] >=0) &&(intensity[i+1][j] <0)) ||((intensity[i][j] >=0) &&(intensity[i][j+1] <0))||((intensity[i][j] <0) &&(intensity[i+1][j] >=0))||((intensity[i][j] <0) &&(intensity[i][j+1] >=0)))
	intensity[i][j] = 255;
      else
	intensity[i][j] = 0;
  for(int i=0;i<maxY;i++)
    intensity[i][0] = intensity[i][1] =  intensity[i][2] = intensity[i][maxX-3] = intensity[i][maxX-2] = intensity[i][maxX-1] = 0;
  for(int j=0;j<maxX;j++)
    intensity[0][j] = intensity[1][j] =  intensity[2][j] = intensity[maxY-3][j] = intensity[maxY-2][j] = intensity[maxY-1][j] = 0;
 
}

void CImage::HoughTransformLine(int nt,int nc,int n1)
{
  Image image(Geometry(nc,nt), "white");
  image.classType(DirectClass);
  Pixels view(image);
  PixelPacket *pixels=view.get(0,0,nc,nt);


  
  int i,j;
  float k;
  int **block = new (int*)[nc+1];
  for(i=0;i<nc+1;i++)
    block[i] = new int[nt+1];
  for(i=0;i<=nc;i++)
    for(j=0;j<=nt;j++)
      block[i][j] = 0;
  float p,maxp;
  maxp = sqrt(maxX*maxX+maxY*maxY);
  float dtheta = 3.14/nt;
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	if(intensity[i][j]>cRes/2)
	  {
	    for(k=-3.14/2;k<3.14/2;k+=dtheta)
	      {
		p = j * sin(k) + i * cos(k);
		p+=maxp;
		block[(int)(p/(2.0*maxp)*nc)][int((k+1.57)/3.14*nt)]++;
	      }
	  }
      }
  bool isline;
  int maxn = 0;
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	if(intensity[i][j]>cRes/2)
	  {
	    isline = false;
	    for(k=-3.14/2;k<3.14/2;k+=dtheta)
	      {
		p = j * sin(k) + i * cos(k);
		p+=maxp;
		if( block[(int)(p/2.0/maxp*nc)][int((k+1.57)/3.14*nt)] >= n1)
		  isline = true;
	      }
	    if(isline)
	      intensity[i][j] = cRes;
	    else
	      intensity[i][j] = 0;
	  }
	else
	  intensity[i][j] = 0;
      }
 for(int i=0;i<nc;i++)
    for(int j=0;j<nt;j++)
      maxn = (maxn>block[i][j]) ? maxn : block[i][j];
  for(int i=0;i<nc;i++)
    for(int j=0;j<nt;j++)
      (*(pixels+i*nt+j))=ColorGray(float(block[i][j])/(float)maxn);
  image.write("temp1.jpg");
  
  for(i=0;i<nc+1;i++)
    delete[] block[i];
  delete[] block;
}

void CImage::HoughTransformCircle(int nc,int n1)
{
  int i,j,k;
  int ***block = new (int**)[nc+1];
  for(i=0;i<nc+1;i++)
    block[i] = new (int*)[nc+1];
  for(i=0;i<nc+1;i++)
    for(j=0;j<nc+1;j++) 
      block[i][j] = new (int)[nc+1];
  
  for(i=0;i<=nc;i++)
    for(j=0;j<=nc;j++)
      for(k=0;k<nc;k++)
	block[i][j][k] = 0;
  float max1,max2,max3;
  max1 = maxX;
  max2 = maxY;
  max3 = sqrt(maxX*maxX+maxY*maxY);

  
  float dx1 = max1/nc;
  float dx2 = max2/nc;
  float i1,j1,k1;
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	if(intensity[i][j]>cRes/2)
	  {
	    for(i1=0;i1<max1;i1+=dx1)
	      for(j1=0;j1<max2;j1+=dx2)
		{
		  k1 = sqrt((i-i1)*(i-i1)+(j-j1)*(j-j1));
		  block[int(i1/dx1)][int(j1/dx2)][int(k1/max3*nc)]++;
		}
	  }
      }
  //  cout<<"done";cout.flush();
  bool iscircle;
  for(i=0;i<maxY;i++)
    for(j=0;j<maxX;j++)
      {
	if(intensity[i][j]>cRes/2)
	  {
	    iscircle = false;
	    for(i1=0;i1<max1;i1+=dx1)
	      for(j1=0;j1<max2;j1+=dx2)
		{
		  k1 = sqrt((i-i1)*(i-i1)+(j-j1)*(j-j1));
		  //		  cout<<int(i1/dx1)<<" "<<int(j1/dx2)<<" "<<int(k1/max3*nc)<<",";cout.flush();
		  if(block[int(i1/dx1)][int(j1/dx2)][int(k1/max3*nc)]>=n1)
		    //	     &&((block[int(i1/dx1)][int(j1/dx2)][int(k1/max3*nc)]>k1*3/2*3.14)))
	
		    {
		      //		      cout<<block[int(i1/dx1)][int(j1/dx2)][int(k1/max3*nc)]<<" "<<k1*3/2*3.14<<" "<<k1<<",";
		      iscircle = true;
		    }
		}
	    if(iscircle)
	      intensity[i][j] = cRes;
	    else
	      intensity[i][j] = 0;
	  }
	else
	  intensity[i][j] = 0;
      }
   cout<<"done";cout.flush();
   /*for(i=0;i<nc+1;i++)
    for(j=0;j<nc+1;j++)
      delete[] block[i][j];
  for(i=0;i<nc+1;i++)
    delete[] block[i];
    delete[] block;*/
  cout<<"end";cout.flush();
}




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);
      }
  return(1);
}

int CImage::saveImage(char* file)
{
  cout<<"0";cout.flush();
  Image image(Geometry(maxX,maxY), "white");
  cout<<"1";cout.flush();
  image.classType(DirectClass);
   cout<<"2";cout.flush();
 Pixels view(image);
   cout<<"3";cout.flush();
 PixelPacket *pixels=view.get(0,0,maxX,maxY);
   cout<<"4";cout.flush();
 int i,j;
  cout<<"5";cout.flush();
  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);
  return(1);
}

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

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

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.  Apply LoG\n";
      cout<<"5.  Find lines - Hough Transform\n";
      cout<<"6.  Find circles - Hough Transform\n";
      cout<<"0.  Exit\n\n";
      cout<<"Enter Option : ";
      cin>>option;
      int div,n1;
      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.LoG();
	  break;
	case 5:
	  cout<<"Enter number of points that must be present in a line : ";
	  cin>>n1;
	  cout<<"Enter the number of divisions in  theta-rough plane : ";
	  cin>>div;
	  i.HoughTransformLine(div,div,n1);
	  break;
	case 6:
	  cout<<"Enter number of points that must be present in a circle : ";
	  cin>>n1;
	  cout<<"Enter the number of divisions in  theta-rough plane : ";
	  cin>>div;
	  i.HoughTransformCircle(div,n1);
	  break;
	case 10:
	  i.print();
	  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.";
	}
    }
}