/***********************************************************
 *       Program Assignment 1---Problem 2
 *          Convolution with Gaussian Mask---Smoothing
 *
 *      Programming with use of the CVIP's Image DATA STRUCTURE: Image
 *
 * Usage:   testDisplay  <input image file> <output image file>
 *
 */


#include <math.h>
#include <stdio.h>

#include <CVIPtoolkit.h>
#include <CVIPdef.h>
#include <CVIPimage.h>
#include <CVIPdrawimg.h>
#include <CVIPspfltr.h>
#include <CVIPview.h>
#include <CVIPobject.h>
//#include <CVIPconvert.h>
#include <CVIPfeatures.h>
#include <histogram.h>
#include <ObjectContour.h>
#include "CVIPgeometry.h"
#include "bilinear.h"
#include "threshold.h"
#include "CVIParithlogic.h"
#include "histogram.h"
#include "CVIPcolor.h"
#include "CVIPhisto.h"
#include <CVIPspfltr.h>
#include <float.h>
#include "CVIPfs.h"
#include <limits.h>
#include "CVIPnoise.h"
#include "CVIPenhance.h"

//#include "Gauss_masks.c"

#define VIEWER "picture"

#define pi 3.141592653589793
void Gauss(float, int, float*);
int main(int argc, char **argv)
{
 Image *origImage, *cvipImage, *cvipHistogramImage, *tempI, *tempI1;
 IMAGE_FORMAT format;
 unsigned iRows, iCols;

 int i, j, k;

 //for Gaussian mask
 float gMask[5];
 float sigma = 1.0;
 int   gMaskSize = 5;
 Matrix *pMatrix;     //pointing to image's pixel matrix
 int sum = 0;
 
 //get the Gaussian parameters
 printf("-----Input Gaussian sigma: "); scanf("%f",&sigma);
 printf("-----Input Gaussian mask size:  "); scanf("%d", &gMaskSize);

 //read in image
 setDisplay_Image(VIEWER, "default");
 format = getFormat_CVIP(argv[1]);
 origImage = read_Image(argv[1], 1); 
 cvipImage = duplicate_Image(origImage);

 //without losing generality, image to be processed 
 //be in unsigned type
 remap_Image(cvipImage, CVIP_BYTE, 0, 255);
 
 //get image size
 iCols = getNoOfCols_Image(cvipImage);
 iRows = getNoOfRows_Image(cvipImage);

 //set image matrix(pixel value) pointer
 //suppose the image be "gray" type
 pMatrix = (Matrix *) cvipImage->image_ptr[0];

 //get the Gaussian mask
 Gauss(sigma, gMaskSize, gMask);
 for (i = 0; i < gMaskSize; i++)
   printf(" %f\n", gMask[i]);


 //Convolve the image with Gaussian mask
 //FIRSTLY, do with Every ROW
 for (i = gMaskSize / 2; i < iRows - gMaskSize / 2; i++)
   for (j = gMaskSize /2; j < iCols - gMaskSize / 2; j++) {
     for (k = 0; k < gMaskSize; k++) {
       sum += pMatrix->rptr[i][j - gMaskSize / 2 + k] * gMask[k];
     }
     sum /= gMaskSize;
     pMatrix->rptr[i][j] = sum; 
     sum = 0;
   }
 sum = 0;
 //SECONDLY, do with Every COLUMN
 for (j = gMaskSize / 2; j < iCols - gMaskSize / 2; j++)
   for (i = gMaskSize /2; i < iRows - gMaskSize / 2; i++) {
     for (k = 0; k < gMaskSize; k++) {
       sum += pMatrix->rptr[i - gMaskSize / 2 + k][j] * gMask[k];
     }
     sum /= gMaskSize;
     pMatrix->rptr[i][j] = sum; 
     sum = 0;
   }
 
 //writing the convolved image
 remap_Image(cvipImage, CVIP_BYTE, 0, 255);
 write_Image(cvipImage, argv[2],CVIP_NO,CVIP_NO,format,1);

 return 0;
} 


//produce Gaussian Mask in H, with INPUTs :sigma = s
//                                         mask size = Hsize
void Gauss(s,Hsize,H)
int Hsize;
float s,*H;
{
  int     i;
  float  cst,tssq,x,sum;
 
  cst=1./(s*sqrt(2.0*pi)) ;
  tssq=1./(2*s*s) ;
 
  for(i=0; i<Hsize; i++) {
    x=(float)(i-Hsize/2);
    H[i]=(cst*exp(-(x*x*tssq))) ;
  }
 
  sum=0.0;
  for(i=0;i<Hsize;i++)
    sum += H[i];
  for(i=0;i<Hsize;i++)
    H[i] /= sum;
}