#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <cv.h>
#include <highgui.h>

//
// general data structures
//

template<class T>
class Image
{
   private:
   IplImage* imgp;
   public:
   Image(IplImage* img=0) {imgp=img;}
   ~Image(){imgp=0;}
   void operator=(IplImage* img) {imgp=img;}
   inline T* operator[](const int rowIndx) {
     return ((T *)(imgp->imageData + rowIndx*imgp->widthStep));}
};

typedef struct{
   unsigned char b,g,r;
} RgbPixel;

typedef struct{
   float b,g,r;
} RgbPixelFloat;

typedef struct{
   float x,y;
} GradVec;

typedef Image<RgbPixel>       RgbImage;
typedef Image<unsigned char>  BwImage;
typedef Image<RgbPixelFloat>  RgbImageFloat;
typedef Image<GradVec>        GradImage;


//
// global variables
//

const int N=4;     // the number of pyramid levels 
int refreshFlag=1; // indicate whether the display needs to be refreshed




void gaussianSmoothing(IplImage *img,float sigma)
{
   RgbImage imagA(img);

   //
   // WRITE CODE FOR SMOOTHING THE IMAGE USING TWO CONVOLUTIONS WITH A 1D GAUSSIAN
   //
}


void computeGradientVectors(IplImage *img, IplImage *grad, IplImage *gradMag)
{
   BwImage   imgA(img);
   GradImage gradA(grad);
   BwImage   gradMagA(gradMag);
   int i,j;

   //
   // REPLACE THE FOLLOWING DEMO CODE WITH YOUR CODE FOR COMPUTING THE
   // GRADIENT VECTORS. COMPUTE THE PARTIAL DERIVATIVES USING A CONVOLUTION
   // WITH A 1D GAUSSIAN IN ONE DIRECTION AND A CONVOLUTION WITH THE
   // DERIVATIVE OF A 1D GAUSSIAN IN THE SECOND DIRECTION
   //

   float theta;
   for(i=0;i<img->height;i++) for(j=0;j<img->width;j++){
     theta=atan2(0.5*img->height-i, j-0.5*img->width);
     gradA[i][j].x=cos(theta);
     gradA[i][j].y= -sin(theta);
     gradMagA[i][j]=int(i*j/(float)(img->height*img->width)*255);
   } 
}


void selectGradientsForDisplay(IplImage *grad, IplImage *gradMag, float prcnt)
{
   GradImage gradA(grad);
   BwImage   gradMagA(gradMag);
   int i,j;

   //
   // REPLACE THE FOLLOWING DEMO CODE WITH YOUR CODE FOR SELECTING X% OF THE
   // GRADIENT VECTORS WITH THE HIGHEST MAGNITUDE
   //

   int step=int(100.0/(prcnt+1));
   if(step<10) step=10;
   for(i=0;i<grad->height;i++) for(j=0;j<grad->width;j++)
     if(i%step || j%step) gradA[i][j].x=gradA[i][j].y=0.0;
}


void computeImagePyramid(IplImage *cimg, IplImage* ip[N], float sigma)
{
   RgbImage cimgA(cimg);
   BwImage ipA[N];
   IplImage *tmpImg;
   BwImage tmpImgA;
   int i,j,k,ik,jk;

   for(k=0;k<N;k++){

     // compute the desired image size
     ik=(int)(cimg->height/pow(2,(N-1)-k));
     jk=(int)(cimg->width /pow(2,(N-1)-k));
     printf("Allocating pyramid level: %d x %d\n",ik,jk);

     // allocate the image
     ip[k]=cvCreateImage(cvSize(jk,ik),IPL_DEPTH_8U,1);
     ipA[k]=ip[k];
   }

   // create the lowest level in the pyramid (grayscale)
   for(i=0;i<ip[N-1]->height;i++) for(j=0;j<ip[N-1]->width;j++)
     ipA[N-1][i][j]=(uchar)(cimgA[i][j].b*0.114 + cimgA[i][j].g*0.587 + cimgA[i][j].r*0.299);

   // initialize the image pyramid
   for(k=N-2;k>=0;k--){

     // clone the previous level and smooth it
     tmpImg=cvCloneImage(ip[k+1]);
     tmpImgA=tmpImg;
     gaussianSmoothing(tmpImg,sigma);

     // sample the smoothed clone
     for(i=0;i<ip[k]->height;i++) for(j=0;j<ip[k]->width;j++)
       ipA[k][i][j]=tmpImgA[i*2][j*2];

     // release the clone
     cvReleaseImage(&tmpImg);
   }

}


void drawGradientVectors(IplImage *grad, IplImage *img)
{
   GradImage gradA(grad);
   RgbImage  imgA(img);
   int lineWidth=1;
   int lineScale=8;
   CvPoint s,e;
   int i,j,scaleDiff;

   scaleDiff=int(img->height/grad->height);

   // draw a green lines to represent the vectors
   for(i=0;i<grad->height;i++) for(j=0;j<grad->width;j++)
     if(gradA[i][j].x!=0 && gradA[i][j].y!=0){

     s.x=j*scaleDiff+scaleDiff/2;
     s.y=i*scaleDiff+scaleDiff/2;
     e.x=int(s.x+gradA[i][j].x*lineScale);
     e.y=int(s.y+gradA[i][j].y*lineScale);

     cvLine(img,s,e,cvScalar(0,255,0),lineWidth);
   }
}


void zoomBwImage(IplImage *src, IplImage *dst)
{
   int i,j,iFac,jFac,ii,jj;
   BwImage srcA(src);
   BwImage dstA(dst);

   iFac=dst->height/src->height;
   jFac=dst->width/src->width;

   // expand src as necessary to fit inside dst
   for(i=0;i<src->height;i++) for(j=0;j<src->width;j++)
     for(ii=0;ii<iFac;ii++) for(jj=0;jj<jFac;jj++)
     dstA[iFac*i+ii][jFac*j+jj]=srcA[i][j];
}

void trackbarHandlerPL(int pos)
{
   refreshFlag=1; 
}

void trackbarHandlerPP(int pos)
{
   refreshFlag=1; 
}

void help(void)
{
   printf("\nUsage: ass1 [<image-file>]\n\n");
   printf("Key summary:\n\n");
   printf("<ESC>: quit\n");
   printf("h    : show this help\n");
   printf("s    : save the current image into a file\n");
   printf("1    : show the color image\n");
   printf("2    : show the grayscale image\n");
   printf("3    : show the gradient magnitude (normalized)\n");
   printf("4    : show the gradient vectors\n");
   printf("\n");
}


int main(int argc, char *argv[])
{
   int displayMode=1;   // the current display mode
   int pyramidLevel;    // the current pyramid level
   int pixelPercentage; // the percentage of pixels in which the
                        // image gradient should be displayed
   float sigma=1.0;     // the standard deviation of the Gaussian filter
   IplImage* cimg = 0;  // the original color image
   RgbImage  cimgA;     // access wrapper for the original color image
   IplImage* ip[N];     // the image pyramid
   BwImage  ipA[N];     // access wrapper for the image pyramid
   IplImage* grad;      // the gradient map
   IplImage* gradMag;   // the gradient magnitude
   IplImage* tmpImg;    // temporary image
   IplImage* zoomImg;   // zoomed image
   int i,j,k,key;

   // capture an image from a camera or read it from a file
   if(argc<2){
     CvCapture* capture = cvCaptureFromCAM(0);
     if(cvGrabFrame(capture)) cimg=cvRetrieveFrame(capture);
   }
   else cimg=cvLoadImage(argv[1]);

   // check the read image
   if(!cimg){
     printf("Could not read/grab image\n\7");
     exit(0);
   }
   cimgA=cimg;
   zoomImg=cvCreateImage(cvSize(cimg->width,cimg->height),IPL_DEPTH_8U,1);

   // create a window with two trackbars
   cvNamedWindow("win1", CV_WINDOW_AUTOSIZE);
   cvMoveWindow("win1", 100, 100);
   pyramidLevel=3;
   cvCreateTrackbar("pyramid level", "win1", &pyramidLevel ,3, trackbarHandlerPL);
   pixelPercentage=10;
   cvCreateTrackbar("pixel percentage", "win1", &pixelPercentage ,20, trackbarHandlerPP);

   // create the image pyramid
   computeImagePyramid(cimg,ip,sigma);


   // enter the keyboard event loop
   while(1){

     key=cvWaitKey(10); // wait 10 ms for a key
     if(key==27) break;

     switch(key){
     case '1':
     case '2':
     case '3':
     case '4':
       displayMode=key-'0';
       refreshFlag=1;
       break;
     case 's':
       if(!cvSaveImage("out.jpg",ip[pyramidLevel])) printf("file save error\n");
       break;
     case 'h':
       help();
       break;
     }

     // update the display as necessary
     if(refreshFlag){
       refreshFlag=0;
       switch(displayMode){

       case 1:
  cvShowImage("win1",cimg);
  break;

       case 2:
         zoomBwImage(ip[pyramidLevel], zoomImg);
  cvShowImage("win1",zoomImg);
  break;

       case 3:
       case 4:
  // allocate the gradient images
         grad=cvCreateImage(cvSize(ip[pyramidLevel]->width,ip[pyramidLevel]->height),
     IPL_DEPTH_32F,2);
         gradMag=cvCreateImage(cvSize(ip[pyramidLevel]->width,ip[pyramidLevel]->height),
     IPL_DEPTH_8U,1);

  // compute the gradient vectors
  computeGradientVectors(ip[pyramidLevel],grad,gradMag);

  // display the gradient magnitude
  if(displayMode==3){
    zoomBwImage(gradMag, zoomImg);
    cvShowImage("win1",zoomImg);
  }

  // display the gradient vectors
  if(displayMode==4){
    selectGradientsForDisplay(grad,gradMag,pixelPercentage);
    tmpImg=cvCloneImage(cimg);
    drawGradientVectors(grad,tmpImg);
    cvShowImage("win1",tmpImg);
    cvReleaseImage(&tmpImg);
  }

  cvReleaseImage(&grad);
  cvReleaseImage(&gradMag);
  break;
  break;
       }
     }

   }

   // release the images
   cvReleaseImage(&cimg);
   cvReleaseImage(&zoomImg);
   for(k=0;k<N;k++) cvReleaseImage(&ip[k]);
   return 0;
}