/**
 * erica bolan
 * convolve.cpp
 * convolve [theta1] [theta2]
 * performs convolution on a given image. if theta value(s) are input, performs
 * gaussian or difference of gaussian convolution; otherwise defaults to the
 * kernel set in the global variables section. the kernel size may be reset
 * for the gaussian filter, but must remain at 3x3 for all other filters.
 * requires a file called "dump.dat" which consists of a single row (or
 * column) of grayscale color values. writes the pixels for the convolved
 * matrix out to "newimage.dat," in a single row of grayscale pixel values.
 * this image can then be viewed with octave or gimp.
 *
 * octave commands:
 * imshow(reshape(dump,320,240)')
 * imshow(reshape(matrix,320,240)')
 * imshow(reshape(newimage,320,240)')
 *
 * simple kernels from http://www.roborealm.com/help/Convolution.php
 */

#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <math.h>
#include <vector>
#include <cstdlib>
using namespace std;

/**global constants*/
const int WIDTH=320;         //original image dimensions
const int HEIGHT=240;
const int k=3;               //size of kernel (square)

/**global variables*/
int image[HEIGHT][WIDTH];    //original image

double theta1,theta2;
double kernel2[k*k];         //only filled for DoG filters
//double kernel1[k*k]={0,-1,0,-1,5,-1,0,-1,0};         //sharpen
//double kernel1[k*k]={-1,-1,-1,-1,8,-1,-1,-1,-1};   //detect lines
double kernel1[k*k]={-1,0,0,0,1,0,0,0,0};          //enhance edges
//double kernel1[k*k]={-2,-1,0,-1,1,1,0,1,2};        //emboss
//double kernel1[k*k]={1,-2,1,-2,4,-2,1,-2,1};       //find edges


/**
 * Writes a pixel matrix to the given file.
 * @param fname name of file to write to
 * @param vals pixels to write (written in column-major format)
 */
void printToFile(char *fname,int vals[HEIGHT][WIDTH])//vals[WIDTH][HEIGHT])
{
  FILE *fp;
  fp=fopen(fname,"w");
  for(int i=0; i<HEIGHT; i++)
    for(int j=0; j<WIDTH; j++)
      fprintf(fp,"%d ",vals[i][j]);
  fclose(fp);
}

/**
 * Fills the kernel with gaussian distribution values.
 * @param theta the value of theta to use in gaussian calculation
 * @param which 1=fill kernel1, 2=fill kernel2
 */
void findGaussian(double theta,int which)
{
  double G, secondTerm, firstTerm;
  double e=2.718281828;
  int x,y,count=0;
  for(x=-k/2; x<=k/2; x++)
  {
    for(y=-k/2; y<=k/2; y++)
    {
      secondTerm=pow(e,-(x*x+y*y)/(2*theta*theta));
      firstTerm=1/(2*3.14159*theta*theta);
      G=secondTerm*firstTerm;
      if(which==1)
	kernel1[count]=G;
      else
	kernel2[count]=G;
      count++;
    }
  }
}

/**
 * Multiplies a given block of the image by the convolution filter.
 * @param pix the current block of pixels to be convolved
 * @param filter the filter (kernel) to be used
 */
int blockMultiply(int pix[k][k], double filter[k*k])
{
  double sum=0;
  for(int i=0; i<k; i++)
    for(int j=0; j<k; j++)
      sum+=pix[i][j]*filter[i*k+j];
  return (int)sum;
}

/**
 * Calculates the matrix for the convolved image.
 * @param which 0=regular convolution, 1=difference of gaussian
 */
void convolveImage(int which) 
{
  int block[k][k];
  int newImage[HEIGHT][WIDTH];
  for(int i=0; i<HEIGHT; i++)
    for(int j=0; j<WIDTH; j++)
      newImage[i][j]=0;
  int avg1,avg2,newPixel;

  for(int picRow=0; picRow<=HEIGHT-k; picRow++)  //big picture row
  { 
    for(int picCol=0; picCol<=WIDTH-k; picCol++)  //big picture column
    {
      for(int blockRow=0; blockRow<k; blockRow++) //kernel row
	for(int blockCol=0; blockCol<k; blockCol++) //kernel column
	  block[blockRow][blockCol]=image[picRow+blockRow][picCol+blockCol];
      avg1=blockMultiply(block,kernel1);
      //difference of gaussian--does the subtraction cause the value to cross 0?
      if(which==1)
      {
	avg2=blockMultiply(block,kernel2);
	newPixel=avg1-avg2;
	if((avg1<0 && newPixel>0) || (avg1>0 && newPixel<0))
	  newPixel=255;
	else
	  newPixel=0;
      }
      //otherwise just make sure they're between 0 and 255
      else
      {
	if(avg1<0)   
	  avg1=0;
	if(avg1>255)
	  avg1=255;
	newPixel=avg1;
      }
      //set the new pixel
      newImage[picRow+k/2][picCol+k/2]=newPixel;
    }
  }

  printToFile("newimage.dat",newImage);
}

/**
 * Reads in the file and stores the values into a matrix on which
 * convolution can be performed. The values are assumed to be grayscale,
 * 0-255. Values must be in a single line of text, with spaces between.
 */
void readFile()
{
  ifstream fin;
  fin.open("dump.dat");
  if(fin.fail())
  {
    cerr << "can't find dump.dat" << endl;
    exit(1);
  }
  int x,row=0,col=0;
  int count=0;

  //FILE *fp;
  //  fp=fopen("vals.dat","w");
  for(int i=0; i<HEIGHT; i++)
  {
    for(int j=0; j<WIDTH; j++)
    {
      fin>>image[i][j];
      //  fprintf(fp,"%d ",image[i][j]);
    }
  }
  fin.close();
  //fclose(fp);


  printToFile("matrix.dat",image);
}

int main(int argc, char *argv[])
{
  readFile();
  if(argc==1)  //simple convolution; default to kernel defined in global variables
  {
    cout<<"simple convolution filter, k="<<k<<endl;
    convolveImage(0); 
  }
  else if(argc==2)   //single gaussian filter; take theta from command line
  {
    theta1=atof(argv[1]);
    cout<<"gaussian filter, k="<<k<<", theta="<<theta1<<endl;
    findGaussian(theta1,1);
    convolveImage(0);
  }
  else if(argc==3)   //difference of gaussian filter; take both thetas from command line
  {
    theta1=atof(argv[1]);
    theta2=atof(argv[2]);
    cout<<"dog filter, k="<<k<<", theta1="<<theta1<<", theta2="<<theta2<<endl;
    findGaussian(theta1,1);
    findGaussian(theta2,2);
    convolveImage(1);
  }
  return 0;
}